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Absttapt 

This report describes a simple and fast atmospheric correction algorithm used to 
correct radiances of scattered sunlight measured by aircraft and/or satellite above a uniform 
surface. The atmospheric effect, the basic equations, a description of the computational 
procedure, and a sensitivity study are discussed. The program is designed to take the 
measured radiances, view and illumination directions, and the aerosol and gaseous 
absorption optical thicknesses to compute the radiance just above the surface, the iiradiance 
on the surface, and surface reflectance. Alternatively, the program will compute the upward 
radiance at a specific altitude for a given surface reflectance, view and illumination 
directions, and aerosol and gaseous absorption optical thicknesses. The algorithm can be 
applied for any view and illumination directions and any wavelength in the range 0.48 |im - 
2.2 pm. The relation between the measured radiance and surface reflectance, which is 
expressed as a function of atmospheric properties and measurement geometry, is computed 
using a radiative transfer routine. The results of the computations are tabulated in a look-up 
table which forms the basis of the correction algorithm. The algorithm can be used for 
atmospheric corrections in the presence of a rural aerosol. The sensitivity of the derived 
surface reflectance to uncertainties in the model and input data is discussed. 
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LQ Introduction 


1.1 Background 

The purpose of this report is to describe a fast and simple atmospheric correction 
algorithm to derive the surface reflectance, and other parameters, from radiances measured 
by satellite and/or aircraft in the visible and near IR parts of the spectrum. The original 
version of this algorithm was developed to correct the radiances measured during FIFE 
(Eirst ISLSCP Eield Experiment - Sellers et al., 1988) which is taking place at the Konza 
Prairie in Kansas. While the algorithm has been designed to correct radiances measured over 
a rural site for the wavelength range 0.48 Jim £ X <. 2.2 pm, a sensitivity study has 
shown that the algorithm can be a practical tool for many applications of remote sensing, for 
which a uniform surface can be assumed, and for which the optical characteristics of the 
aerosol do not differ significantly from the rural aerosol (the algorithm should not be applied 
for correcting for the effect of desert dust or fog). As a result we would like to bring the 
algorithm to the attention of the scientific and engineering community. 

1.2 The Intervening Atmosphere 

Atmospheric aerosols, which are liquid or solid particles suspended in the air, have a 
significant importance in evaluating satellite imagery for remote sensing of the earth's 
surface. The atmospheric aerosol results from natural sources (e.g. desert dust, 
condensation and oxidation of gases released from the biosphere and oceans) and 
anthropogenic sources (e.g. biomass burning and the industrial emission of gases which 
participate in atmospheric chemical reactions and condense into liquid particles). In the 
Southern Hemisphere, near Australia, the aerosol concentration is usually very low (aerosol 
optical thickness in the visible is less than 0. 10) due to the low population, large ocean 
areas, and low humidity. In desert areas, dust storms can increase the optical thickness (t) to 
x = 2.0 and above (hiding the sun). In the Northern Hemisphere the concentration may be 
rather large during long period of times, (e.g. the aerosol optical thickness is around 0.6 in 
the Eastern part of the United States during July and August - Kaufman and Fraser, 1983; 
Peterson et al., 1981), due to industrial pollution. In the state of Rondonia, Brazil, biomass 
burning due to deforestation generates dense smoke (optical thickness 1. 0-3.0) that covers 
the area during most of the dry season. For an aerosol optical thickness over land larger than 
0.20, aerosols affect a significant share of the outgoing visible radiation for a cloudless sky. 
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Therefore, since aerosols affect satellite imagery of the earth's surface, attempts for its 
correction should be taken. 

Satellite images of the Earth's surface in the solar spectrum are contaminated by 
sunlight scattered towards the sensor by atmospheric molecules, aerosols, and clouds (path 
radiance). In addition, solar energy that is reflected from the Earth's surface and serves as 
the remote sensing signal, is attenuated by the atmosphere. This combined atmospheric 
effect is wavelength dependent, varies in time and space, and depends on the surface 
reflectance and its spatial variation. Correction for this atmospheric effect can produce 
remote sensing signals that are more closely related to the surface characteristics. Molecular 
scattering and absorption in the atmosphere can be accounted for satisfactorily. Gaseous 
absorption is minimized by choosing sensor bands in atmospheric windows. Therefore, 
aerosol scattering and absorption, and the presence of subpixel clouds, are the main 
variables in the atmospheric effect on satellite imagery. 

For a cloudless sky, aerosol scattering is the major variable component of the 
atmospheric effect for dark surfaces, while aerosol absorption is important for bright 
surfaces (Fraser and Kaufman, 1985). In order to perform atmospheric corrections of 
remotely sensed data, the optical characteristics of the atmosphere must be estimated. These 
characteristics may be given in varying levels of detail-from considerable detail (profile of 
the extinction coefficient, the single-scattering albedo and the scattering phase function), to 
less detail (the vertical optical thickness, the average aerosol scattering phase function and 
single-scattering albedo). 

In order to demonstrate the effect of the atmosphere on remote sensing we shall 
discuss vegetation as an example. We may distinguish among the following ecological 
regions: 

- Over remote land areas, with no substantial anthropogenic aerosol contribution, and 
no dust, the aerosol optical thickness may vary between 0.02-0.10. For this variation and 
for typical aerosol characteristics (single scattering albedo G)o = 0.96), the reflectance of the 
atmosphere alone will increase about 0.01 (Fraser and Kaufman, 1985), and a vegetation 
index (ratio between the difference in reflectances between the near IR and the visible and 
the sum - NDVI) of NDVI=0.60, as measured by satellite, will change to 0.58 (Holben, 
1986). These changes are relatively small and, therefore, atmospheric correction in this case 
is not important for most applications. (Remote sensing of ocean color is affected by even a 
small aerosol optical thickness (Gordon et al., 1983).) 
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- Over typical land areas, anthropogenic aerosols and/or dust may generate optical 
thicknesses in the range 0.05-0.25. The corresponding atmospheric effect would change a 
surface reflectance of p = 0.02 to 0.04, and vegetation index of NDVI = 0.60 to 0.55. 
These are significant errors which necessitate atmospheric correction. 

- Over polluted areas, with anthropogenic aerosols from industrial sources (e.g all of 
Eastern U.S. and Europe during the summer) or areas affected by dust, fog or smoke 
(tropical regions, regions in the far east and Sahel), the aerosol optical thickness may vary in 
the range 0.1 -1.0. In this range of variation the atmospheric effects are very large. The 
surface reflectance would vary from 0.02 to 0.08 and the vegetation index would decrease 
from 0.6 to 0.45. 

1.3 Atmosp heric Corrections 

The correction procedure requires information about the atmospheric optical 
characteristics. Due to the difficulty in determining these characteristics, the only operational 
use of atmospheric corrections today is that of the ocean color (Gordon et al., 1983), where 
the corrections depend on the condition of the very low reflectance of the water in the red. 
Otherwise, information on the atmospheric optical characteristics can be obtained from three 
different sources: 

Climatology: Documented information on the atmospheric characteristics and their 
variation can be used to estimate the expected atmospheric effect for a specific part of the 
world and a specific season. Such documentation can be obtained from the analysis of 
measurements taken from the ground, and partially from the analysis of satellite data (Fraser 
et al., 1984; Kaufman, 1987; Kaufman et al., 1988). This source of information will be 
used for optical characteristics that cannot be determined otherwise for the particular image 
being corrected. 

Measurements from the ground: The aerosol optical thickness can be obtained 
from sun-photometer measurements (King et al., 1978; Kaufman and Fraser, 1983). The 
phase function can be determined from inversion of solar almucantar measurements, and the 
single- scattering albedo from the collection of particles on filters, preferably by aircraft 
sampling of the entire atmospheric boundary layer. The single-scattering albedo can also be 
determined by measurements of the diffuse and direct flux (Herman et al., 1975; King and 
Herman, 1979; King, 1979) and by lidar techniques (Spinhime et al., 1980). The 
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application of such measurements for atmospheric corrections is useful for intense field 
measurements, or for establishing the climatology of a given area. 

Determination from satellite imagery: For the purpose of atmospheric 
corrections, the path radiance and the corresponding aerosol optical thickness can be derived 
from radiances detected by the satellite over a dark surface. Examples include many land 
surfaces in the blue spectrum, dense dark vegetation in the visible channels (Kaufman and 
Sendra, 1988), and water in the red and near IR. The wavelength dependence of the derived 
aerosol optical thickness (when available) can be used to estimate the particle size and the 
scattering phase function. In the past satellite imagery has been used to determine the 
aerosol optical thickness and other aerosol characteristics. The aerosol optical thickness has 
been derived from satellite imagery of oceans (Griggs, 1975; Mekler et al., 1977; Carlson, 
1979; Koepke and Quenzel, 1979; Takayama and Takashima, 1986), and recently, over 
dense dark vegetation (Kaufman and Sendra, 1988). By using the difference in the 
brightness between a clear and a hazy day, Fraser et al. (1984) demonstrated that the 
difference in the optical thickness can be derived where the surface reflectance is less than 
0.1. Determination of the aerosol single-scattering albedo and particle size was suggested by 
Kaufman et al. (1988) and Kaufman (1987), and applied to trace the evolution of smoke 
from a large forest fire (Ferrare et al., 1988). This method is useful to determine the aerosol 
characteristics (from imagery that includes water-land interfaces) in areas that suffer from 
substantial aerosol outbreaks (e.g. desert dust storms, smoke from fires and concentrated 
anthropogenic aerosol). 

Once the atmospheric characteristics are specified, an atmospheric correction can be 
performed with an equation relating relating measured radiance to the optical properties of 
the atmosphere and surface. The computation requires application of complex and time 
consuming radiative transfer programs (Dave, 1972 a,b,c,d; Ahmad and Fraser, 1982). This 
report presents an algorithm that simplifies the correction procedure by using an a priori 
prepared look-up table that is based on radiative transfer computations. In essence, the 
algorithm simplifies the atmospheric correction procedure to a desktop operation, by 
sacrificing the flexibility to select specific aerosol size distributions and refractive indexes, 
but not optical thickness. 
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1.4 The algorithm 


The algorithm is designed to compute the upward radiance for a given surface 
reflectance or to compute the surface reflectance for a given measured radiance, for almost 
any wavelength in the visible and near-IR spectrum (with appropriately specified gaseous 
absorption), for a wide range of observation zenith angles, solar illumination angles and 
azimuth angles between the observer and the solar rays, as well as any height of the 
observer (aircraft or satellite). Any practical value of the aerosol optical thickness can be 
used, but the algorithm is restricted to a specific aerosol size distribution and refractive 
index. 

The relation between the measured radiance and the surface reflectance is expressed as a 
function of the path radiance, downward flux at the ground, atmospheric transmission, and 
the atmospheric backscattering ratio. Using this relation, a look-up table is constructed 
which relates the measured upward radiance to surface reflectance for several aerosol optical 
thicknesses, solar zenith angles, measurement wavelengths, and a range of observation 
directions. This look-up table is based on the tabulation of the results of radiative transfer 
computations which are made using a Dave (1972 a,b,c,d) code. It is assumed that the 
atmosphere and surface are horizontally homogeneous, and the surface reflects light 
according to Lambert's law. The light scattered by the atmosphere and the surface is 
assumed to be unpolarized. The atmosphere is also assumed to be cloud-free. 

Radiation properties of the cloudless atmosphere depend on both molecular and aerosol 
constituents. Molecular scattering and absorption, except for water vapor absorption, are 
easy to account for. Aerosol effects are more variable and are therefore more difficult to 
correct. The parameters used to describe these aerosol effects are: the optical thickness 
which determines the amount of extinction, the single scattering albedo which determines the 
fraction of light scattered from the total extinction, and the single scattering phase function 
which describes how the light is scattered as a function of direction. For the most part, 
aerosol extinction is the dominant parameter in the aerosol component of the atmospheric 
effect (Fraser and Kaufman, 1985). Thus, in this algorithm the aerosol optical thickness is 
the only variable aerosol parameter. The algorithm uses a constant aerosol single scattering 
phase function and scattering albedo chosen to represent a rural environment. Because these 
assumptions can introduce error in the derived surface reflectances, a sensitivity analysis is 
performed to estimate the uncertainty in in derived surface reflectances. 
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In the first part of this document, the atmospheric effect and its effects on the surface 
reflectance are discussed. Next, a description of the equations used in the algorithm and the 
construction of the look-up table is given. A sensitivity study is then performed to estimate 
the errors associated with the initial assumptions, interpolation algorithm, and uncertainties 
in the input data. Finally, the FORTRAN code for the algorithm is listed. 

2.0 Atmospheric Effect 

2.1 General Discussion 

The atmospheric effect is caused by the scattering and absorption of solar radiation 
by molecules and aerosols. There are three components to this effect : 

1) The downward solar radiation is absorbed and scattered by the atmosphere and it is 
diffused by forward scattering. Because this diffusion increases the angular distribution of 
the radiation, the downward radiation interacts with the surface in a wide range of 
directions. Thus the surface reflection coefficient for this radiation is different from the 
reflection coefficient for the direct solar beam (Lee and Kaufman, 1986). 

2) Radiant energy reaching a remote sensor is reflected from the ground both within the 
instantaneous-field-of-view (ifov) and from the region outside of it. Part of the reflected 
energy within the ifov is transmitted directly to a sensor and can be considered signal; the 
remaining radiation is absorbed and scattered. Part of the radiation that is reflected from 
outside of the ifov passes through the column containing the ifov and is scattered there 
towards the sensor. This component is associated with the adjacency effect. It augments 
the measured radiance and is partially corrected for when deriving the surface reflectance. 

3) Radiation is scattered by the atmosphere into the ifov without being reflected by the 
surface. This component is called the path radiance and increases the apparent reflection of 
the surface. 

An example of the difference between the spectral surface reflectance for typical 
vegetation (assumed Lambertian) and the corresponding upward radiance above the 
atmosphere is shown in Figure 1. In this figure, the upward radiance is normalized by the 
incident solar flux to produce reflectance units. This normalized radiance is the apparent 
reflectance as seen from the sensor. The reflectance of the earth-atmosphere system can be 
greater than the surface reflectance, the same, or weaker. In the visible spectrum the 
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• RADIANCE /F 0 /i. 0 



WAVELENGTH (/im) 


Figure 1. Surface reflectance ( ) and the corresponding radiance ( — ) above the 

atmosphere for a typical vegetation. The radiance is normalized by Fq, the incident solar 
flux, and by Po» the cosine of the solar zenith angle (from Kaufman, 1987). 
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surface reflectance is weak, the path radiance is relatively strong, and the reflectance above 
the atmosphere exceeds that at the ground. In the near infrared (0.7<X<1.6 pm) the 
surface reflectance is strong, and the reflectances at the surface and above the atmosphere 
are nearly the same. The loss of radiation from the surface by extinction is augmented at 
the same rate by atmospheric scattering. For wavelengths longer than 1.6 |im, atmospheric 
scattering does not compensate for the attenuation loss. 

2.2 Mathematical Description 

In order to correct for the atmospheric effects discussed in the previous section, a 
relation is developed between the upward spectral radiance L m measured from satellite or 
aircraft and the surface reflectance p: L m = f(p). The radiance is equivalent to the specific 
intensity as defined by Chandrasekhar (1960, p.l), except that the radiance, as used here, 
is the radiant energy within a unit wavelength interval instead of the energy per unit 
frequency. The function f depends on the atmospheric and surface optical properties, 
observation and sun directions, and wavelength. The radiance L m can be expressed 
explicitly as a function of the path radiance L 0 (upward radiance for zero surface 
reflectance), the downward flux through a horizontal surface at the ground (for zero 
surface reflectance), the total (direct + diffuse) transmission from the surface to the 
observer T, and the atmospheric backscattering ratio s. It is assumed that the atmosphere 
and the surface are horizontally homogeneous, but the atmospheric optical properties vary 
in the vertical direction. The surface is assumed to reflect light according to Lambert's law. 
The light scattered by the cloud-free atmosphere and surface is assumed to be unpolarized. 
The relation between L m and p is (Chandrasekhar, 1960) 


L m = L + 
o 


(PFd'fl 

K (1 - S p) 


( 1 ) 


Here L m is the spectral radiance measured from aircraft or satellite and is a function of X* 
e o> x a> x gs> w o> *g. 0 , Z, and <p , where 


X is the wavelength of the radiation, 

0 Q is the solar zenith angle, 

t is the aerosol optical thickness (used with base e), 
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x gs is the molecular scattering optical thickness, 

to o is the ratio of the aerosol scattering and extinction optical thicknesses, 

x g is the gaseous absorption optical thickness, 

Z is the observation height, 

0 is the propagation direction zenith angle of the radiant energy at the ground, 

cp is the azimuth angle (azimuthal angles are measured with respect to the 

principal plane through the sun; 0° lies in the plane containing the direction 
of propagation of die direct sunlight). 


Figure 2 shows the angular coordinates used in the algorithm. The functions L 0 , T, F d and 
s have the following functional dependances: 


L 0 = L 0 (X,e o ,x a ,x gs ,x g ,co 0 ,Z,e,(p) 
T = T(X,x a ,x gs ,x g ,co o ,Z,0) 


F d = F d^’ 0 O’ X a’ X g S ’V COo) 

S = s(A.,X a ,X gs ,tg,C0o) 


The correction algorithm is based on the inverse of eq. (1), where the surface reflectance p 
can be expressed in terms of the measured radiance L m : 


f 

(1+sf) 


( 2 ) 


where 


*(L m -L o ) 


(3) 


The algorithm will compute L m using (1) if p is given or will compute p using (2) and (3) 
if L m is given. Other quantities computed are: the total irradiant flux F g on a horizontal 
surface at the ground, 


F 

g 


(1 -sp) 


(4) 
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Zenith 

4 


Sun 



Figure 2. Angular coordinates used in the algorithm. The X-Y plane is a horizontal plane 
tangent to the earth's surface at the observation point. The solar zenith angle 0 O , 
observation zenith angle 0, observation scan angle 0’, and observation azimuth angle <J> are 
shown. In this particular representation of planar geometry 0 = 0’; the general relationship 
is given by eq. (9). 


10 




and the upward radiance L g at the ground in the direction of observation. 


L 

g 


p p d 

7C (1 - S p) 


( 5 ) 


3.0 Correction Algorithm 

The correction program is baspd on the tabulation of the results of 
radiative transfer computations of 1^, F d , s, and T/n (note that the 
transmission T is not tabulated). The primed parameters are normalized flux 
and radiance rather than their absolute values. The normalized values are 
related to their corresponding absolute values by the following equations: 


F . (0 , t ) = 
d v o’ a 7 


F cos 0 
o o 


, kL 

l„ <»„. 0. 4>. V Z) = V <«) 

F cos 0 
o o 


The value of F 0 represents the solar spectral flux passing through a surface orthogonal to 
its propagation at the top of the atmosphere. In order to use the table for atmospheric 
correction, the measured absolute radiance L m is converted to L m ' using 



t m 

n L 

F o' cos 0 o 


(7) 


where 


F ' 
o 




( 8 ) 


and d is earth-sun distance for the day of the year when measurements are made, d is the 
mean earth-sun distance, and F 0 is the solar flux now computed for each of the spectral 
bands using the solar spectral flux data from Neckel and Labs (1984). Since the look-up 
table tabulates observation zenith angle 0 of the line-of-sight from the ground to the 
observing platform, the scan angle 0' measured by aircraft or satellite is converted to 
observation zenith angle 0 using 


11 



( 9 ) 


0 = sin' 1 (l+f)sin0’ 
r s 


where Z is the height of the sensor above the ground, and r s is the radius of the earth. 

The computations are performed by a Dave code (1972a, b, c, d) with a series of 
radiative transfer programs. These programs compute the flux and radiance of the scattered 
radiation emerging at any level of a plane-parallel atmosphere. Henceforth, primed fluxes 
and radiances indicate that they have been normalized as in eq. 6. Variables with the 
superscript m are measured. 


3.1 Model Wavelengths 


The look-up table is computed for the following wavelengths: 0.639, 0.845, 0.486, 
0.587, 0.663, 0.837, 1.663, and 2.189 pm which correspond to the following sensors: 
NOAA-9 AVHRR band 1 (0.58 - 0.68 pm), band 2 (0.725 - 1.10 pm); Landsat-5 TM 
and NS-001 TMS band 1 (0.45 - 0.52 pm), band 2 (0.52 - 0.60 pm), band 3 (0.63 - 
0.69 pm), band 4 (0.76 - 0.90 pm), band 5 (1.55 - 1.80 pm), and band 7 (2.10 - 
2.35 pm). These wavelengths are listed in Table 1. The wavelength chosen to represent a 
particular band is computed by first calculating X* 


, f X L m F T'dX. 
X = J- ■ 7 ? 

j L m F q 'F dX 

where 


( 10 ) 


L m ' = normalized radiance at the top of the atmosphere 
F 0 = extraterrestrial solar spectral flux 


'F = response function of the sensor 

The values for F 0 are obtained from Neckel and Labs (1984) while the values for 'F 
correspond to the specific sensor. For the NOAA-9 AVHRR and Landsat TM sensors, 
these values are obtained from Kidwell (1985) and Markham and Barker (1985), 
respectively . The effective wavelength A,* in eq. 10 is a function of the normalized 
radiance, which is a function of the surface reflectance, aerosol optical thickness, and 
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Table 1. Spectral bands, aerosol refractive indices, and optical thicknesses. 



AVHRR 



TM 




1 

2 

1 

2 

3 

4 

5 

7 

Sensor Wavelengths 









50% minimum response (nm) 

569.8 

714.3 

452.4 

528.0 626.4 776.4 1567.5 2097.2 

50% maximum response (nm) 

699.3 

982.2 

517.8 

609.3 693.2 904.5 1784.1 2349.0 

peak response (nm) 

680.0 

760.0 

503.0 

594.0 677.0 800.0 1710.0 2200.0 

model (equation 7) (nm) 

639.0 

844.6 

486.2 

586.9 662.7 837.3 1662.7 2188.6 

Indices of Refraction 









n' (accumulation mode) 

1.43 

1.43 

1.43 

1.43 

1.43 

1.43 

1.40 

1.40 

k (accumulation mode) 

io- 8 

io- 8 

io- 8 

10 8 

10 8 

io- 8 

10 4 

10 4 

n’ (coarse particle mode) 

1.53 

1.53 

1.53 

1.53 

1.53 

1.53 

1.40 

1.35 

k (coarse particle mode) 

10-7 

10-7 

10-7 

10-7 

10-7 

10-7 

10- 4 

0.00814 

Optical Thicknesses 









Molecular Scattering x e . 

0.0540 

1 0.0180 

0.159 

0.0841 

0.0449 

0.0176 

0.0012 

0.0004 

Gaseous Absorption 









Ozone Tg ®3 

0.0240 

0.00064 0.00663 0.0317 

0.0174 

0.0000 

0.0000 

0.0000 

Water Vapor Tg H 2° 

0.00605 0.0933 

0.00000 0.0002 

0.0068 

0.0410 

0.0957 

0.0741 

CO 

Carbon Dioxide x g 2 

0.00071 0.0146 

0.0000 

0.0000 

0.0000 

0.0021 

0.0077 

0.0091 

Aerosol Absorption (x a = 0.25) 








xa a 

0.0148 

0.0223 

0.0130 

0.0142 

0.0150 

0.0218 

0.0423 

0.0189 

Composite (section 3.31 









T H 

X g 

0.0247 

0.0152 

0.0066 

0.0317 

0.0174 

0.00206 0.00771 0.00908 

T L 

0.0208 

0.1156 

0.0130 

0.0142 

0.0218 

0.0628 

0.138 

0.0930 


0.0455 

0.131 

0.0196 

0.0459 

0.0392 

0.0648 

0.146 

0.102 
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geometry. The normalized radiance is computed for different models of the earth- 
atmosphere system representative of surface and atmospheric conditions expected for the 
FIFE Konza Prairie site. Figure 3 shows the surface reflectance profiles used to model this 
site. A similar approach can be used to determine the wavelengths which correspond to 
other sensors. 

The algorithm is generalized to accept any wavelength in the range 0.48-2.2 pm. For 
wavelengths for which there are no entries in the look-up table (see Table 1), the algorithm 
will interpolate the atmospheric functions (Lq, T, F^, s) for the desired wavelength. The 
interpolation is performed assuming that aerosol parameters are proportional to the 
wavelength raised to a power: 

In Tg S _ - In X, In x a ~ - In X, In P _ In X, 

where P is the scattering phase function. As a result, the radiances in the look-up table can 
be interpolated linearly between the wavelengths on a log - log scale. The only nonlinear 
relation is between the gaseous absorptions in the different wavelengths. Correction for the 
gaseous absorption is discussed in section 3.3. 

3.2 Aerosol Properties 

Because the absorption and scattering of light by atmospheric aerosols is highly 
variable, some assumptions regarding the size, shape, and composition of the aerosol must 
be made. In this model, the aerosols are assumed to be spheres so that Mie theory can be 
used to calculate the scattering by aerosols. Although in general the aerosol particles are not 
spherical, it is assumed that the sizes assigned to the aerosol particles are the sizes of 
spheres that have similar scattering properties to the measured aerosol distribution (Shettle 
and Fenn, 1979). This assumption has further basis because it has been found that the 
aerosol particles become more spherical as the relative humidity increases (Nilsson, 1979). 

The algorithm uses a bimodal aerosol size distribution which combines the optically 
effective fraction of the accumulation mode (0. 1 pm < d <; 1.0 pm) and the coarse 
particle mode (d > 0.5 pm). For the accumulation mode, the dry particles are assumed to 
be composed of 80% water soluble sulfates and 20% water insoluble, dust-like material 
(Nilsson, 1979). At 70% relative humidity, water composes half of the volume of these 
aerosols. The coarse particle mode aerosols are assumed to be made of mostly water 
insoluble, dust-like particles. The size distribution of the aerosols is represented as the sum 
of two log-normal distributions; the two distributions represent the accumulation and coarse 
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Figure 3. Surface reflectance for senescent grass from a burned surface on the Konza 

Praire at the FIFE site. The first profile ( ) is for nadir view, measured at 14:38 CDT 

on 15 July 1986 while the second profile ( — ) is for an observation zenith angle of 45° 
measured at 15:13 CDT on 15 July 1986 (Asrar, 1986, private communication). 
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particle modes. The number density function of particles of radius r per cubic centimeter of 
air per micrometer of radius is (Shettle and Fenn, 1979) 


n (r) 


= dN(iO = Y — 
dr In 

i=l L 


N. 

i 


- (log r - log rj,) 2 ~| 

i In (10) r a. -/Tic 

exp 

to 

-a, 

i 


( 11 ) 


where 


N(r) = cumulative number density of particles of radius r 
Oi = standard deviation of the logarithm of the radius 
= geometric mean radius 
Nj = total number density in i* mode 

The values of Ni, r n ‘, and Ci used in the model correspond to the 70% relative humidity, 
rural aerosol model of Shettle and Fenn (1979); these values are shown in Table 2. The 
values of Ni shown in Table 2 are normalized such that Ni + N2 = 1 particle/cm 3 . This 
bimodal aerosol size distribution is assumed constant with height. 

The composition of the aerosols is expressed in terms of the complex refractive index 
n = n' - i k. The refractive index for both the accumulation and coarse particle modes is 
assumed to depend on the wavelength. Five different real refractive indices are chosen 
(Nilsson, 1979) and are listed in Table 1. The imaginary index of refraction is modeled 
assuming that most of the aerosol consists of weakly absorbing particles with k < 1(H> 
mixed with a small number of highly absorbing particles with k ~ 1.0. Although Shettle 
and Fenn (1979) choose to represent mixtures of this type with a composite imaginary 
refractive index in the range 0.001 < k < 0.01, which is close to the values obtained by 
various remote sensing and in-situ techniques (Patterson and Grams, 1984; Reagan et al., 
1980), this procedure is not adopted here because as Bohren and Huffman (1983) point 
out, no common substances exist which have an imaginary index in this range. In the 
correction algorithm, the imaginary refractive index corresponding to the weakly absorbing 
particles is used. The imaginary refractive indices for the accumulation mode, which are 
listed in Table 1 , correspond to water while the values for the coarse particle mode 
correspond to crystalline quartz, which is a constituent of atmospheric dust (Nilsson, 
1979). The aerosol refractive index is assumed to be constant with height. 
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Table 2. Aerosol size distribution parameters. 


Accumulation mode Coarse particle mode 


Geometric mean radius r n * 

0.0285 

0.457 

Standard deviation of the 



logarithm (base 10) 

0.81 

0.81 

of the radius aj 



Number density Nj 

0.999875 

0.000125 
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The aerosol absorption needed to match the experimentally measured value is modeled 
by adding the necessary absorption to the gaseous absorption values. Thus, it is assumed 
that the absorbing particles are small compared to the wavelength and that they occur 
separately from the other particles (external mode); in this case their scattering effects are 
small relative to their absorptive effects (Fraser and Kaufman, 1985). The required 
additional aerosol absorption amount is determined by the aerosol single scattering albedo 
C0o which is the ratio of scattering to extinction. Aerosol absorption is then given by 1 - coq. 
The values of (Do used in the algorithm are derived from the 70% relative humidity, rural 
aerosol model of Shettle and Fenn (1979). The single scattering albedos, shown in 
Table 1, range from 0.95 at 0.486 |xm to 0.86 at 2.550 |im. These values for the visible 
spectrum agree with those reported by Waggoner et al. (1981) which lie in the range 
0.89 £ cd 0 £ 1.0 for rural areas. The amount of additional aerosol absorption which is 
added to the gaseous absorption is 

x a = x (to ' - 0 ) ) (12) 

a a v o o' v ' 


where 


t a = aerosol optical thickness 

CDo' = single scattering albedo for chosen refractive index 

(Do = desired single scattering albedo 

Since in general k « 1, then cd 0 ' = 1, and 

\ - x (1 -cd ) (13) 

a a o 

Should the algorithm be applied for an urban aerosol or smoke, where the single scattering 
albedo cd q used in the formation of the look-up tables (see Table 1) may be smaller or larger 
than the new value oo 0 *, the difference in absorption can be corrected by adding (or 
subtracting) the excess absorption to the gaseous absorption: 

At^ = T a (CD q - CD*) (14) 


3.3 Gaseous Absorption 

For the most part, the Landsat TM and NOAA AVHRR visible and near IR bands 
have been selected to minimize gaseous absorption. However, in some cases, the sensor 
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channel is either relatively broad (as in the case of AVHRR band 2) or lies within a broad 
gaseous absorption band (as in the case of TM bands 2 and 3 which lie within the ozone 
continuum) so that the absorption by atmospheric gases can be both significant and 
variable. In the atmospheric correction algorithm, gaseous absorption was computed using 
the LOWTRAN 6 code (Kniezys et. al., 1983) which computes atmospheric absoiption 
from 0.250 |im to 28.5 Jim due to water vapor, carbon dioxide (and other uniformly 
mixed gases), ozone, nitrogen continuum, oxygen, and HNO 3 . In the case of the AVHRR 
and TM bands, the absorbing gases are water vapor, carbon dioxide, and ozone. The 
gaseous absorption optical thickness in the vertical direction due to gaseous species x for 
each band is computed using 



J T x L m F 0 'P dX 

jL^Fo'RdX 


(15) 


where 

m = air mass along inclined path 

T x = transmittance due to gaseous species x 

Since the gaseous absorption is quite variable, the algorithm uses a weighted average of the 
absorption values computed using the tropical, mid-latitude summer, and mid-latitude 
winter atmospheres given in the LOWTRAN 6 code. The weighted average of the gaseous 
absorption x* used in the algorithm is given by 

5 

T x * = 0.25 x\ + 0.5 T x + 0.25 x X (16) 

g gl g2 g3 

XX X 

where Xgj, x^ and x^ are the gaseous absoiption values computed using the mid-latitude 
winter, mid-latitude summer, and tropical values respectively. The absorption optical 
thicknesses due to water vapor, carbon dioxide, and ozone for each band are shown in 
Table 1. 

The algorithm can be applied to any specified gaseous optical thickness Xg. An 
approximate correction is applied to the radiances in the look-up table to account for the 
excess (or deficit) in the absorption (AXg). The user may compute the required value of Xg 
based on the sensor spectral response using the LOWTRAN program (Kneizys et al., 

1983). For the correction of the look-up table we may distinguish between AXgL - excess 
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or deficit in the gaseous absorption in the lower part of the atmosphere (e.g. water vapor), 
assumed to be mixed uniformly with the aerosol; and Aigfj excess or deficit in the gaseous 
absorption in the upper part of the atmosphere, above the aerosol layer (e.g. O 2 , 03 , 

CO 2 ). The correction is performed by the following transformations in the look-up table: 


L o“* L o ex P 


Ax - 


(17a) 


F d^ F d ex P 


-< At gL + At gH>' 




(17b) 


T — >T exp 


< A V + At gH> 


(17c) 


s— >sexp[ - 2 Ax g J 


(17d) 


Here p = cos 0 and Po = cos 0 O . In these equations the effect of multiple scattering on 
the path length through the atmosphere is neglected. It is also assumed that the path 
radiance L 0 is generated above the middle of the boundary layer. As result, the additional 
gaseous attenuation is made by half of the boundary layer, but all of the atmosphere above. 
The parameter s is the reflectance of the atmosphere for radiation entering its base. The 
effective direction for reflection is 60°. Hence, the effective absorption optical thickness is 
twice the vertical value. Theses assumptions are based on our physical understanding of 
radiative transfer and are tested in Section 6 . 

3.4 Altitude Profiles 

The radiative transfer computations require as input the altitude distributions of both 
the absorbing gases and aerosols. These two distributions have separate profiles. The 
altitude distribution of aerosols used in the algorithm is based on the 'average' distribution 
of Braslau and Dave (1973) which is shown in Figure 4. The aerosol altitude distribution is 
first scaled to obtain the desired aerosol optical thickness. 
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Figure 4. Aerosol particle number density profiles for: 1) 'average' distribution of Braslau 
and Dave (1973) used in the correction algorithm; 2) 50 km surface visibility , 3) 23 km 
surface visibility, and 4) 10 km surface visibility models from Shettle and Fenn (1979). 
Figure 4a is for altitudes between 10 and 70 km while figure 4b is for altitudes between 0 
and 10 km. These profiles have been normalized for an aerosol optical thickness x a = 0.25 
at 0.55 iim. 
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Because the various gaseous absorbers described in the previous section have 
different altitude distributions, a composite altitude distribution is computed which accounts 
for all the absorbing gases including the aerosol absorption described in section 2.2. In this 
method, the total gaseous absorption is divided into "low" and "high" components. The 
"low" component is composed of water vapor and aerosol absorption while the "high" 
component is composed of ozone and CO 2 absorption. The "low" component uses an 
altitude distribution based on the 'average' aerosol profile of Braslau and Dave (1973), 
while the "high" component uses an altitude distribution based on the mid-latitude ozone 
profile of McClatchey et al. (1971). The altitude distributions of the "low" and "high" 
components are normalized and combined to produce an altitude profile that retains the 
maximum aerosol and water vapor absorption near the surface and the maximum ozone 
absorption in the lower stratosphere. Because the aerosol absorption is a function of the 
aerosol optical thickness (see eq. 13), the gaseous absorption profile is a function of 
aerosol optical thickness as well as wavelength. The atmospheric pressure profile used in 
the algorithm is adapted from the mid-latitude summer profile of McClatchey et al. (1971). 

3.5 Molecular and Aeros ol Optical Thickness 

The molecular scattering (or Rayleigh) optical thickness t gs is computed for sea 
level from (Hansen and Travis, 1974) 

T gs = 0.008569 X 4 (1 + 0.01 13 X 2 + 0.00013 X ' 4 ) (18) 

where the wavelength is in micrometers. This expression assumes the surface sea-level 
pressure is 1013 mb. In the case the surface is not at sea level. The value of the molecular 
scattering optical thickness t gs is assumed to vary according to: 

v (Z 0 > = tgs (°) exp (-^) (19) 

where Zq is the height of the surface above sea level in kilometers. In order to avoid a need 
for a new radiative transfer computation for each height of a surface, the computations are 
performed for Zq = 0.4 km, and the look-up table is adjusted each time the user specifies a 
different height The algorithm adjusts the look-up table to account for the different 
molecular optical thickness by adjusting the wavelength of the radiation. Substitution of the 
relationship between X and 'tg (eq. 18) yields 
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„ (Zo - 0.4) 

X(Zq) = \(0.4) exp f 36 “ ] 


( 20 ) 


For example, for Zo = 0 km, X will increase by 1.1%, whereas for Zq = 2 km, X will 
increase by 4.3%. The error introduced by this method is in the effective scattering phase 
function, but this change is negligible relative to the general uncertainty in the scattering 
phase function due to the uncertainty in the aerosol size distribution. 

The relation between upward radiance and surface reflectance for each spectral band 
is computed for four aerosol optical thicknesses, x a = 0.0, 0.25, 0.50, and 1.00, except for 
TM bands 5 and 7 where only the first two values are used. These values are selected to 
cover the range of aerosol optical thicknesses which could be expected for most remote 
sensing applications. 

3.6 Measurement Altitudes 

The radiative transfer computations are tabulated at three measurement altitudes 
above the ground: 0.45 km, 4.5 km , and 80 km. The algorithm will use these altitudes to 
interpolate to the input measurement altitude. Since 80.0 kilometers is above the 
atmosphere, corrections to satellite measurements are made with the 80 km tables. 

For altitudes less than 0.45 km, it is assumed that the aerosols and the absorbing gases 
are well mixed and it is possible to linearly interpolate between the atmospheric optical 
properties (Lq and T) at Z = 0 km and at Z = 0.45 km above the ground: 

Lo (Z) = L 0 (Z = 0.45) (q^) (21a) 

T (Z) = (T(0.45) Z + (0.45-Z)) / 0.45 (2 1 b) 

For altitudes above 4.5 km, interpolations are performed between values of L 0 and 
T at 4.5 km and 80 km. Although the atmospheric model contains aerosols between these 
altitudes, the interpolations are based on the assumption that Lo and T depend essentially on 
molecular scattering between 4.5 and 80 km. The interpolations are performed linearly as a 
function of the molecular scattering coefficient (jgs(Z): 

o gs (Z) = o gs (0) exp(^) (22) 

Hence the transmission and Lq become 
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T(Z) = (T(80.0) - T(4.5)) 


exp(-4.5/9) - exp(-Z/9) 
exp(-4.5/9) - exp(80/9) 


(23) 


L 0 (Z) = (L o (80.0) - L 0 (4.5)) 


exp(-4.5/9) - exp(-Z/9) 
exp(-4.5/9) - exp(80/9) 


(24) 


Between 0.45 km and 4.5 km the interpolation is uncertain due to the inhomogeneity of 
the aerosol layer. Since the linear interpolation used below 0.45 km is usually suitable to 
extrapolate for heights below 1 km, and the aerosol concentration decreases rapidly above 
3 km so that the interpolation used above 4.5 km is appropriate for extrapolation to heights 
between 3 and 4.5 km, the algorithm extrapolates from these two regions to the desired 
height h (0.45 km < Z < 4.5 km), and chooses the result that shows a minimal atmospheric 
effect. Thus, the lower values of L 0 and higher values of T are chosen for the solution. 

4.0 Look-up Tables 

The look-up tables used by the correction algorithm contain the normalized radiances 
and fluxes computed by the radiative transfer routines. The values generated are stored to 
be used by the main program. Eight look-up tables are produced, one table for each of the 2 
AVHRR bands of 0.639 and 0.845 pm and one for each of the 6 TM bands of 0.486, 
0.587, 0.663, 0.837, 1.663 and 2.189 pm. Each table is arranged for 3 heights of 0.45, 

4.5 and 80.0 kilometers. Each table contains values for 9 solar zenith angles 0 O (10, 20, 
30, 40, 50, 60, 66, 72, and 78°), 13 observation zenith angles 0 (0° to 78°, every 6°), 19 
observation azimuth angles (0° to 180°, every 10°, plus 5° and 175°), and 4 aerosol optical 
thicknesses t a (0.0, 0.25, 0.50 and 1.0) for all wavelengths, except for 1.663 and 
2.189 pm where only first two optical thicknesses of 0.0 and 0.25 are used. 

5.0 Computational Procedure 

5.1 General 

Input data to the correction program consists of the date, time, measurement 
wavelength, aerosol and gaseous absorption optical thicknesses at the measurement 
wavelength A., solar zenith angle 0 O , observation scan angle 0', observation azimuth angle 
<)>, height of the surface above sea level Zq, measurement height Z m , and the measured 
spectral radiance in absolute L m or reflectance L m ' units (option 1), or surface reflectance p 
(option 2). The wavelength, solar and observation angles, and aerosol optical thickness 
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data must be in the ranges discussed above as no extrapolation is performed. If the selected 
wavelength or altitude does not match the values used to construct the look-up table, the 
algorithm interpolates on wavelength and altitude as described above. 

If option 1 is selected, the program computes the surface reflectance p, total spectral 
irradiance on the surface F , and total spectral radiance of the ground L in the direction of 

© 5 

observation using eq. 2, 3, 4 and 5. The total spectral irradiance F is computed in 
Watts/m 2 /pm and total upward spectral radiance L is computed in Watts/m 2 /|im/sr. If 

o 

option 2 is selected, the program computes the absolute and normalized radiances L m and 
L m ', total spectral irradiance on surface Fg and the total upward spectral radiance L g in 
Watts/m 2 /|im/sr. 

The program for making atmospheric corrections consists of a main program called 
FIFEWAV. The computations are performed in single precision and seven subroutines 
are called at different stages by the main program. The listings of the main program and 
subroutines are given in section 9. The statement numbers for the program and subroutines 
are given in parentheses. Table 3 lists the variables appearing in the main program and 
associated subroutines and their equivalent in the text. The asterisk superscript indicates 
interpolated values. The subroutines called by the main program FIFEWAV are listed 
below: 

Subroutine READIN (6350 - 6580): This subroutine reads the input data from unit 
number 5 and writes it on unit number 6. The purpose of this subroutine is to perform a 
check on the input data set. 

Subroutine FINDW (6620 - 6860) : This subroutine checks to see if falls in the 
range of wavelengths for which look-up tables are available and picks two wavelengths 
between which X m falls. The subroutine prints an error message and returns the control to 
the main program if A, m does not fall in the range of wavelengths look-up table provides. 
The main program processes a new data point. 

Subroutine INTSFX (6900 - 7100): This subroutine computes values of the 
sun-earth distance R for 365 days of the year. These values are returned to the main 
program, which chooses the value of R for the day the measurements were made. 

Subroutine INTHGH (7140 - 8530): This subroutine computes inteipolated values 
of Lo and T/it for the measured height MHGHT (eq. 21a , 21b, 23 and 24). The 
transmission T is divided by it to account for the it appearing in eq. (1). 
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Table 3. Variables appearing in FIFEWAV program and 
associated subroutines, and their equivalent in the text. 


VARIABLES 

IN 

PROGRAM 

ABSRAD 

AIRR 

AMUO 

ANGLE, PHI 
ARAD 

FDOWN 

FFLUX 

FINT 

FOIRR 

FT 

INT 

IRRID 

MHGHT 

MINT 

MPHI 

MTAU 

MTHET 

MTHETO 

MWAV 

NFDOWN 

OPTH 

PIT 

R, RR 
RHO 

RHOS 

SBAR 

THE 

THETO 

WAV 


VARIABLES 

IN 

TEXT 


Lg 

(section 5) 

V 

(section 5) 

Cosine (0 o m ) 

(eq. 1, 22) 

<P 

(eq. 1) 

V 

(section 5) 

* 

f d 

(section 5) 

i 

F o 

(eq. 22, 23) 

1 * 

L o 

(section 5) 

F o 

(eq. 6) 

T*/jc 

(section 5) 

L'o 

(section 5) 


(section 5) 

z m 

(eq.l) 

L m > L m 

(eq. 1, 22) 

cp m 

(eq. 1) 

^a m 

(eq. 1) 

0 m 

(eq. 1) 

a m 

e o 

(eq. 1) 

X m 

(eq. 1) 

■ * 

F d 

(section 5) 


(eq. 1) 

T/tc 

(section 5) 

R 

(eq. 23) 

P 

(section 5) 

* 

P 

(section 5) 

s 

(eq. 1) 

0 

(eq. 1) 

e 0 

(eq. 1) 

X 

(eq ■ 1) 


26 



Subroutine INTERP (8570 - 9020): The purpose of this subroutine is to return to the 

'4* Jfc 

main program the interpolated values of L c , F d , T /tc, and s*. This is a general 
purpose interpolation routine and is called at different stages by the main program. This 
subroutine does not allow any extrapolation and an error message is printed when the 
values are out of bounds. The asterisk indicates in most cases values interpolated from the 
look-up table. 

Subroutine INTEXP (9040 - 9840): This subroutine sends exponantially interpolated 

*4« *4c % * 

values of L Q , F c j, T /tc and s for the measured optical thickness MTAU. It calles a 
systems subroutine ZXGSN residing in IMSL math package. 

Subroutine ZXGSN (9520): This subroutine is called by subroutine INTEXP, it 
computes the minimum of function on certain interval which is used for exponantial 
interpolation. 

PROCEDURE: 

Subroutine READIN is called (490) to read the input data from unit number 5 and 
write it on unit number 6 as well as to perform a check on the input data. The main program 
reads the first two data cards from unit number 5 and writes the labels for output on unit 
number 6 and 56 (690 - 840). The program reads from unit 20 the solar spectral 
irradiances, which are used to compute the solar flux for the measured wavelength MWAV 
(870 - 950). The program enters into a loop over the number (INUM) of data points to be 
processed (1020). The program reads the values of time, X m , x a m , 0 o m » 0 m , <p m , Z m , L m 
or L m ' , XgL, XgH, and Zq if NOPT is 1, and values of time, X m , x a m , 0 o m * 0 m , cp m , Z m , 
P m > XgL, XgH and Zo if NOPT is 2 (1 190 - 1230). If values of absorptions XgL and XgH 
are 0.0, default values are computed by interpolating linearly from the values of 
TAUGL(X) and TAUGH(X) ( 330 - 340) for the measured wavelength X m (1300 - 1600). 

Subroutine FINDW is called to select the look-up table for the two wavelengths 
between which X m lies by assigning the proper file number NFILE1 and NFILE2 (1700 
- 1740). A new data point is read if subroutine FINDW does not find the input 
wavelength lying in the range of wavelengths for which look-up table is provided (0.486 - 
2.2 [im). 

The ratio of earth-sun distance (R, eq. 8) is computed for 365 days of the year by 
subroutine INTSFX. The program (1870 - 2020) changes the month and day to the 
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Julian day if MOPT is 1 . The corrected values of incident solar flux F Q are computed by 
using eq. 8 (2150). 

The observation scan angle of satellite or aircraft is changed to observation zenith angle 
at the ground by using eq. 9 (2250). 

Statements (2610 - 3040) read the look-up table chosen for the two wavelengths on 

it 

each side of the measured wave length X™ 1 . The parameters 0, 9 , X, x & , 0 O , s, F d , L Q , 
and T/jc are read for four optical thicknesses, except for wavelengths of 1.663 pm and 
2.189 Jim, for which look-up tables are available for only two optical thicknesses of 0.0 
and 0.25. The subroutine INTHGH (3120) is called to compute the interpolated values of 
Lo and T /it for the measured height MHGHT (eq. 21a, 21b, 23 and 24). The statements 
(3170 - 3180) adjust the look-up table for the surface height SHGHT and then interpolate 
values of s, Fd , L 0 ’ and T /71 (3230 - 3410) for the measured wavelength X m (see 
documentation sec. 3.1). The excess or deficit gaseous absorption in the upper and lower 
atmosphere is computed (3450 - 3500) using TgL and XgH whose values are supplied as 
input parameters; or if values read are zero, default values computed earlier in program 
(1300 - 1600) are used. New values of s, Lo', Fd', and T/tc are computed (3540 - 3780) 
after adjusting for the excess and deficit of gaseous absorption (eq. 17a, 17b, 17c, and 
17d). 

The next step is to compute the interpolated values of Lo', Fd', and T/tc for the 
measured height and geometry (0 o m , 0 m , and <p m ). The mesh of the tables is small 
enough to allow accurate linear interpolations. The interpolations are made for each of the 
four aerosol optical thicknesses: 0.00, 0.25, 0.50, and 1.00. L'* 0 is interpolated from the 
data set Lo' for 0 o m , <p m , and 0 m ; F*d is interpolated from the data set Fd' for 0 o m ; and 
T */k is interpolated from the data set T/tc for 0 m . To make these interpolations, subroutine 
INTERP (4050 - 4140) is called to compute the interpolated values of both Lo’ and Fd’ for 
0 o m . The part of program between 4270 to 4430 calls subroutine INTERP to compute 
the interpolated values of Lo’ for <p m . Subroutine INTERP is called again to compute the 
interpolated value of Lo' and T/tc for 0 m (4680 - 4770). If the measured 0 Q m , cp m , or 0 m 
are out of the range of data tabulated for 0 O , (p, and 0, subroutine INTERP returns control 
to the main program after printing an appropriate message. It does not allow any 
extrapolation. The program will then process a new data point if the error message is 
printed. 
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Subroutine INTEXP (4910) is called to perform exponential interpolations on four 
radiation parameters L"* 0 , F*d» s, and T * In for measured optical thickness MTAU(4960 - 
5190). Subroutine INTEXP checks the four parameters one by one for linearity and sends 
the control back to the main program if any of the functions are linear. The main program 
then calls the subroutine INTERP to perform linear interpolation. 

Statement 5270 sends control to statement 5550 to compute the spectral radiance, if 
NOPT is 2. Statements 5280 to 5380 use the interpolated values L * 0 , F * d , and T * In 
to compute the surface reflectance p, total spectral irradiance F , and total spectral radiance 

o 

of the ground L in the direction of observation (eq. 2, 3, and 4). The total spectral 

o 

irradiance F g is computed in Watts/m 2 /|im, and the total upward spectral radiance L g is 
computed in Watts/m 2 /|im/sr. Statements 5440 - 5450 write the output on FORTRAN 
logical units 6 and 56. 

Statements 5540 - 5790 use the interpolated values L* 0 , F *d> and T*/n to compute the 
spectral radiance L m , total spectral irradiance F g , and total spectral radiance of the ground 
L in the direction of observation (eq. 2, 3, and 4) if MOPT is 2. The total spectral 

5 

irradiance F g is computed in Watts/m 2 /|im, and the total upward spectral radiance L g is 
computed in Watts/m 2 /pm/sr. The output is written on FORTRAN logical units 6 and 56. 

5.2 Input 

Input data to this program consists of two input files. The first file (FORTRAN 
logical unit number 5) describes input to a particular case, while the second file 
(FORTRAN logical unit number 20) reads solar spectral irradiances (Neckel and Labs, 
1984) as function of wavelength. The program uses 8 look-up tables which are not to be 
changed by the user (FORTRAN logical unit numbers 7 - 14). In order to demonstrate the 
input data, four cases representative of different options are presented. The input cards for 
all four runs are listed in Table 4. 

INPUT FROM UNIT NUMBER 5: (see Table 4) 

1. The first card contains the labels to identify the input parameters in step 2 below. 

2. The second card contains the options for determining the units for the measured 
reflectance or radiance, the format for the date of the measured data, and the option either 
to compute the surface reflectance or radiance. 
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TABLE 4. EXAMPLES OF INPUT DATA. 


(RUN 1) 

IOPT MOPT NOPT 
111 

INUM MON DAY YEAR TIME ZONE 

6 7 28 1987 CDT 

TIME MWAV MTAU MTHETO MTHET MPHI MHGHT MINT TAUL TAUH SHGHT 

1010 0.486 0.250 26.00 0.00 0.00 80.00 00.11409 00.0001 00.0070 00.400 

1011 0.486 0.250 26.00 0.00 0.00 4.50 00.08089 00.0001 00.0070 00.400 

1012 0.486 0.250 26.00 0.00 0.00 0.45 00.05125 00.0001 00.0070 00.400 

1013 0.486 0.250 26.00 18.00 150.00 80.00 00.12345 00.0001 00.0070 00.400 

1014 0.486 0.250 26.00 18.00 150.00 4.50 00.08649 00.0001 00.0070 00.400 

1015 0.486 0.250 26.00 18.00 150.00 0.45 00.05263 00.0001 00.0070 00.400 

(RUN 2) 


IOPT 

MOPT ] 

NOPT 









1 

2 

1 









INUM 

DAY 

YEAR TIME ZONE 








6 

209 1987 CDT 








TIME 

MWAV 

MTAU 

MTHETO 

MTHET 

MPHI 

MHGHT 

MINT 

TAUL 

TAUH 

SHGHT 

1010 

0 .486 

0.250 

26.00 

0.00 

0 . 00 

80.00 

00 .11409 

00.0001 

00 . 0070 

00 . 400 

1011 

0.486 

0.250 

26.00 

0.00 

0.00 

4.50 

00.08089 

00.0001 

00.0070 

00.400 

1012 

0.486 

0.250 

26.00 

0.00 

0.00 

0.45 

00.05125 

00.0001 

00.0070 

00.400 

1013 

0.486 

0.250 

26.00 

18.00 

150.00 

80.00 

00.12345 

00.0001 

00.0070 

00.400 

1014 

0.486 

0.250 

26.00 

18.00 

150.00 

4.50 

00.08649 

00.0001 

00.0070 

00.400 

1015 

0.486 

0.250 

26.00 

18.00 

150.00 

0.45 

00.05263 

00.0001 

00.0070 

00 . 400 


(RUN 3) 

IOPT MOPT NOPT 
2 11 

INUM MON DAY YEAR TIME ZONE 


6 7 28 1987 CDT 


TIME 

MWAV 

MTAU 

MTHETO 

MTHET 

MPHI 

MHGHT 

MINT 

TAUL 

TAUH 

SHGHT 

1010 

0.486 

0.250 

26.00 

0.00 

0.00 

80.00 

63.12500 

00.0001 

00.0070 

00.400 

1011 

0.486 

0.250 

26.00 

0.00 

0.00 

4.50 

44.75500 

00.0001 

00.0070 

00 . 400 

1012 

0.486 

0.250 

26.00 

0.00 

0.00 

0 . 45 

28.35600 

00.0001 

00.0070 

00 . 400 

1013 

0.486 

0.250 

26.00 

18.00 

150.00 

80.00 

68.30400 

00.0001 

00.0070 

00.400 

1014 

0.486 

0.250 

26.00 

18.00 

150.00 

4.50 

47.85900 

00.0001 

00.0070 

00.400 

1015 

0.486 

0.250 

26.00 

18.00 

150.00 

0.45 

29.12000 

00.0001 

00.0070 

00.400 

(RUN 

4) 










IOPT 

MOPT NOPT 

1 o 









L 

INUM 

1 4 

MON DAY YEAR TIME 

ZONE 







6 

7 

28 198^ 

r CDT 








TIME 

MWAV 

MTAU 

MTHETO 

MTHET 

MPHI 

MHGHT 

MRHO 

TAUL 

TAUH 

SHGHT 

1010 

0.486 

0.250 

26.00 

0.00 

0.00 

80.00 

00.05000 

00.0001 

00.0070 

00 . 400 

1011 

0.486 

0.250 

26.00 

0.00 

0.00 

4.50 

00.05000 

00.0001 

00.0070 

00.400 

1012 

0.486 

0.250 

26.00 

0.00 

0.00 

0.45 

00.05000 

00.0001 

00.0070 

00.400 

1013 

0.486 

0.250 

26.00 

18.00 

150.00 

80.00 

00.05000 

00.0001 

00.0070 

00.400 

1014 

0.486 

0.250 

26.00 

18.00 

150.00 

4.50 

00.05000 

00.0001 

00.0070 

00.400 

1015 

0.486 

0.250 

26.00 

18.00 

150.00 

0.45 

00.05000 

00.0001 

00.0070 

00.400 


30 



IOPT , MOPT , NOPT (FORMAT: 315) 


For IOPT = 1 
For IOPT = 2 

For MOPT = 1 
For MOPT = 2 
For NOPT = 1 
For NOPT =2 


measured radiance should be in reflectance units - L m '. 

measured radiance should be in absolute radiance units of 
W atts/m 2 /sr/jim. 

the day, month and year should be provided. 

the day of year (Julian date) and year should be given. 

surface reflectance (RHOS) is computed. 

radiance (MINT) is computed for the input geometry and height; the 
measured surface reflectance (MRHO) is input. The variable name 
MINT is used for programming convenience. MINT is a computed 
and not measured when NOPT = 2. 


Run 1 and run 2 choose the IOPT = 1 option few selecting reflectance units of 
measured radiance and NOPT = 1 option to compute surface reflectance. MOPT (the 
option to select the day of year) is set to 1 and 2, respectively. Run 3 chooses IOPT = 2 
and MOPT = NOPT = 1. Run 4 selects NOPT to be 2 and IOPT and MOPT to be 1. 


3. The third card contains the labels to identify the input parameters in step 4 below. 

4. The fourth card contains the information about the number of data points INUM to be 
processed, the date, and information about the time zone (e.g Central, Pacific, or Mountain 
and Daylight or Standard Time). The information about the time zone is not used for any 
computations. There are two options for this input card depending on the value of MOPT 
(FORMAT 415, 30A1): 

For MOPT = 1: INUM, IMONTH, ID AY, IYEAR, TIME ZONE 

e.g.: (1, 12, 25, 1987, CDT) for 1 data point to be processed for 
December 25, 1987, and the time zone is CDT or Central Daylight Time. 

For MOPT = 2: INUM, IDAY, IYEAR, TIME ZONE 

e.g.: ( 1, 359, 1987, CDT) 
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The fourth row of run 1 reads number of data points, month, day, and year. The fourth 
row of run 2 reads number of data points, the Julian day and year for which the 
measurements are made. The radiance is read in reflectance units for both runs 1 and 2. 

5. The fifth card contains the labels to identify the input parameters in step 6 below. 

6. The sixth card can be repeated for the number of data points INUM. The data points 
should be for the same date. The following names in parentheses are the variable names 
used by the program. The sixth card contains values of time in hours and minutes (without 
punctuation; this input parameter is not used in any computations but is merely for the 
record) (MTIME), measured wavelength X m (MWAV) in micrometers, measured 
aerosol optical thickness x a m (MTAU) for the measured wavelength X m , solar zenith angle 
0 o m (MTHETO) in degrees, observation scan angle 0 m (MTHET) in degrees, 
observation azimuth angle <p m (MPHI) in degrees, observation height Z m (MHGHT) in 
kilometers, and spectral radiance L m or L m (MINT) from satellite or aircraft in reflectance 
or absolute units depending on option (IOPT ) read in the first card, TgL (TAUGL), XgH 
(TAUGH), and Zo (SHGHT) in the following format: 

MTIME, MWAV, MTAU, MTHETO, MTHET, MPHI, MHGHT, MINT, 
TAUGH, TAUGL, SHGHT 

FORMAT (15, lx, F 6.3, F 6.3, 4F7.2,F10.5,3F8.4) 

Card 6 for run 3 reads the measured radiance in absolute units . Card 6 for run 4 when 
NOPT is 2 reads the surface reflectance value MRHO in the following format: 

MTIME, MWAV, MTAU, MTHETO, MTHET, MPHI, MHGHT, MRHO, 
TAUGH, TAUGL, SHGHT 

FORMAT (15, lx, F 6.3, F 6.3, 4F7.2,F10.5,3F8.4) 


INPUT FROM LOOK-UP TABLE (UNIT NUMBERS 7 - 14): 

The input from the look-up table is read after the execution of the statements given 
in parentheses. Subroutine FINDW (1610) searches the data base and then reads only the 
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required table for the two wavelengths between which measured wavelength X m (MWA V) 
lies. The following are obtained: 

1. values of observation zenith angle 0 in degrees (2610, 2650) 

2. values of observation azimuth angle <p in degrees (2610, 2660) 

3. values of wavelength X in micrometers, aerosol optical thickness t a , reflectance of 

» 

atmosphere s, flux incident on surface F d (2830 - 2840) 

4. a blank line 

i 

5. values of atmospheric radiance L Q in reflectance units (2920 - 2950) 

6. values of transmission from surface T/7t (3000). 

5.3 Output 

A successful execution of this FORTRAN program will result in the output 
described below, which is written to FORTRAN logical unit numbers 6 and 56. Four 
cases representative of the different options discussed in section 5.2 are shown in Table 5. 

OUTPUT TO UNIT NUMBER 6: (see Table 5.) 

The innput data set is written as it is read from unit number 5. The message 'END OF 
INPUT DATA SET AS READ FROM UNIT NUMBER 5' marks the end of input data set. 
Hie next line of this output contains the information about the date for which measurements 
are made: 

For MOPT = 1, month, day, year, and time zone will be written (Table 5, run 1). 

For MOPT = 2, Julian day, year, and time zone will be written (Table 5, run 2). 

Each subsequent line gives the measured and derived output for each data point processed. 
The first seven entries on a line provide information about the following input parameters 
as read from logical unit number 5: time (Central, Pacific, or Mountain and Daylight or 
Standard) in hours and minutes, measured wavelength (X m ) in micrometers, measured 
aerosol optical thickness (x a m ), solar zenith angle (0 o m ) in degrees at the time of 
observation, measured observation zenith angle (0 m ) in degrees, measured observation 
azimuth angle (<p m ) in degrees, and measured observation height (Z m ) in kilometers. The 
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TABLE 5. EXAMPLES OF COMPUTED PRODUCTS WITH INPUT FROM TABLE 
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1014 0.486 0.25 | 26.00 18.01 150.00 4.50 46.47 0.0865 ] 1650. 26.97 0.051 

1015 0.486 0.25 1 26.00 18.00 150.00 0.45 28.27 0.0526 I 1650. 26.40 0.050 



(RUN2) 

START OF INPUT DATA AS READ FROM UNIT NUMBER 
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TABLE 5. CONTINUED 
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eighth and ninth entries depend on 10 PT. ForlOPT = 1, the measured radiance L m is 
given in reflectance units and labeled under RELATIVE RADIANCE; the corresponding 
radiance in absolute units is computed by the the program and is labeled under 
RADIANCE. For IOPT = 2, the measured radiance L m is given in absolute units and is 
labeled under RADIANCE; the corresponding radiance in reflectance units is computed by 
the the program and is labeled under RELATIVE RADIANCE. Next three surface 
quantities are computed and given as values for total spectral inradiance F g , total spectral 
radiance L , and surface reflectance p (eq.2,3 and 4). The output in Table 5 for run 4 

o 

computes the eighth and ninth columns, since NOPT = 2 and the surface reflectance is 
input (column 12). 

OUTPUT TO UNIT NUMBER 56: (see Table 6) 


The label for this output shown in table 6 reads ' INTERPOLATED RADIATION 
PARAMETERS. The second line gives information about the date and time zone of the 
input data set. Each subsequent line gives the measured and derived output for each data 
point processed. The first seven entries on a line provide information about the following 
input parameters as read from logical unit number 5: time (Central, Pacific, or Mountain 
and Daylight or Standard) in hours and minutes, measured wavelength (A, m ) in 
micrometers, measured aerosol optical thickness (t a m ), solar zenith angle (0 o m ) in degrees 
at the time of observation, measured observation zenith angle (0 m ) in degrees, measured 
observation azimuth angle ((p m ) in degrees, and measured observation height (Z m ) in 
kilometers. The eighth and ninth entries depend on IOPT and NOPT. For IOPT = 1 
and NOPT = 1, the measured radiance L m is given in reflectance units and labeled under 
RELATIVE RADIANCE; the corresponding radiance in absolute units is computed by the 
the program and is labeled under RADIANCE. ForlOPT = 2 and NOPT = 1, the 
measured radiance L m is given in absolute units and is labeled under RADIANCE; the 
corresponding radiance in reflectance units is computed by the the program and is labeled 
under RELATIVE RADIANCE. IF NOPT = 2 the values of radiance in absolute and 
reflectance units are computed by the program. Next three quantities are computed and 
given as values for L 0 > F d ', s, and transmission T (eq. 6). 


5.4 Summary of Error Messages 

A summary of error messages is given in Table 7. First, examples of input data are 
given followed by the different error messages that would result. For the sixth row and 
second column of the upper table, where the input measured wavelength X m = 0.345, it 
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TABLE 1. EXAMPLES OF ERROR MESSAGES FOR VARIOUS ERRORS IN THE INPUT DATA 
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does not lie between the range (0.486 - 2.2 (im ), for which the look-up table is available 
(section 4.0). The subroutine FINDW prints an error message 'Measured wavelength 
0.345 is out of range. The data point #2 is not processed; Actual range 
(0.486 - 2.2 fim)'. For the following rows (2, 3, 4, and 5) of input data, the measured 
aerosol optical thickness x™, solar zenith angle 0™ observation zenith angle 0 m , and 
observation azimuth angle (p 01 arc out of the ranges of values provided by the look up table 
(section 4.0). The subroutine INTERP prints error messages stating the variables which 
are out of range. 

6.0 Sensitivity Study 

A careful error analysis would be extensive, because so many measurement and 
surface parameters arc involved. Instead an attempt is made to estimate the maximum 
errors that might occur during FIFE in estimating the surface reflectance (p). The errors 
given here arc are not root-mean-square errors. The surface radiance L g error can be 
calculated with the reflectance error and eq. 5. The irradiance at the ground (F g (eq.4) and 
also Fd (eq. 4 and 5)) is not affected appreciably by the perturbations. The main source of 
error usually is caused in estimation of the path radiance (Lo). The absolute surface 
reflectance error is insensitive to strength of the reflectance. 

The perturbations are placed in three categories: model errors, interpolation errors, and 
uncertainties in the input data (input errors). The model parameters studied are: 

1) uncertainty in the aerosol single scattering phase function 

2) uncertainty due to the neglect of polarization in the radiative transfer computations 

3) uncertainty in water vapor and ozone absorption 

4) uncertainty in aerosol absorption 

5) uncertainty in the height distribution of aerosols 

Since the algorithm must interpolate parameters on aerosol optical thickness Xa, solar zenith 
angle 0o, observation zenith angle 0, azimuth angle <|>, wavelength X, and height Z, the 
errors introduced by these interpolations are studied. Finally, because the input data will 
have some uncertainty associated with them, the effects of errors in the input radiance, 
input aerosol optical thickness, and input geometry are studied. The studies described 
above are mostly performed for a wavelength X = 0.486 pm, since the atmospheric effects 
generally are the greatest at the shortest wavelength. Unless stated otherwise, the 
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unperturbed simulation model has an aerosol optical thickness Ta = 0.25, three observation 
altitudes (0.45 km, 4.5 km , 80 km), and two view geometries (Model 1: 0o = 30°, 0 = 
0°, <|> = 0°; Model 2: 0o = 30°, 0 = 18°, <t> = 150°). In order to estimate the largest 
errors, the wavelength, geometry, and aerosol optical thickness are varied in a few cases 
where the errors are larger. 

In the discussion which follows, a description of each type of sensitivity test is 
given in Table 8, followed by a discussion of the results. The first part gives results of 
model errors, the second part of interpolation errors, and the third part of input errors. A 
summary of the sensitivity results is shown in Table 9, which shows the errors in the 
derived surface reflectance for two observation geometries in the case of the model errors. 
For interpolation and input errors, the resulting errors in the derived surface reflectance are 
shown for only one geometry. Table 10 gives additional details concerning the errors and 
radiation parameters. 

6.1 Model Errors 

The atmospheric correction algorithm causes errors because the model does not 
represent exactly the state of the atmosphere at the time of an observation. The simulated 
measured remote radiance or surface reflectance is first computed with a radiative transfer 
model that contains the perturbation. These computations produce a given perturbed 
radiance L m + AL m at the sensor altitude, as well as the other radiation parameters. The 
perturbed radiance L m + AL m serves as an input in the correction algorithm, resulting in a 
derived perturbed surface reflectance p + Ap. 

6.1.1 Aerosol Single Scattering Phase Function 

The scattering effects of aerosols increase with increasing aerosol optical thickness, 
which usually is largest at the shortest wavelengths; but compared with molecular 
scattering, the relative contribution by the aerosols to backward scattering decreases with 
decreasing wavelength. Therefore, an intermediate wavelength of 0.639 jxm (AVHRR 
band 1) is used for the phase function sensitivity study. 

The phase function sensitivity test is performed using two phase functions which 
differ appreciably from the one used in the model. As discussed previously, the aerosol 
single scattering phase function used in the algorithm is computed for two log-normal 
aerosol size distributions combined, which are chosen to represent accumulation and coarse 
particle modes.The first phase function used in this sensitivity analysis is chosen from a 
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Table 8. Case numbers assigned to Sensitivity Tests for the Atmospheric Correction Algorithm 
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Upward radiance <Watts/nr /^lm/sr) at the ground for no error 


Input radiance (reflectance units) at sensor height for no error 
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Water vapor absorption optical thickness Tg = 0.0486 (AVHRR band 2) 


vapor absorption optical thickness Tg = 0.1235 (AVHRR band 2) 
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Input error in aerosol optical thickness of Ala - 0.10 for true Ta 



Table 8 continued 
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Table 9. Summary of sensitivity results. 



Table 9 continued 


1 

i 
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0.486 0.1122 -0.1 0.05 0.0 0.1230 -0.1 0.0 -0.0 aerosol height distribution (23 km visibility) 

0.486 0.1122 -0.2 0.05 0.0 0.1230 -0.3 0.05 0.0 aerosol height distribution (10 km visibility) 
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Table 10. Sensitivity Results 
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B) Interpolation errors. 
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Table 10 continued 
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library of phase functions that were computed for 19 models of the accumulation mode. 

The size distribution is a single log-normal aerosol size distribution with an index of 
refraction of n = 1.349 - 0.008i, geometric mean mass radius r m = 0.6 |im, and standard 
deviation of the logarithm of the radius a = 0.61. This phase function is compared with 
the phase function used in the algorithm in Figure 5 along. The second phase function used 
in the sensitivity analysis is chosen assuming a power law aerosol size distribution. The 
power chosen is v = 4 where 


dn 

din r 


~r v 


and the index of refraction is n = 1.53 - 10* 7 i. This phase function is also shown in 
Figure 5. Chayanova and Shifrin (1966) found this power law most closely matched the 
model number 4 (the surface visibility is 20 km) aerosol phase function measured by 
Barteneva (1960). These two models are designated as cases 1 and 2 in Tables 8-10. 

Figure 5 shows that the bimodal log-normal phase function used in the algorithm lies 
between the two phase functions used in the sensitivity study for scattering angles greater 
than about 50°. Therefore, these two phase functions are chosen to represent extreme 
experimental conditions where the true aerosol scattering in the backward direction does not 
match the scattering assumed in the model. 

Both of the phase functions used in the sensitivity analysis produced errors in the 
derived surface reflectance of approximately ±20-30%. The phase function derived using 
the single mode log-normal aerosol size distribution produced surface reflectance values too 
small, while the phase function derived using the power size distribution produced surface 
reflectance values too large. It should be noted that the two phase functions were chosen to 
represent the possible extremes in the scattering phase function (for large scattering angles) 
so that the errors associated with uncertainties in the phase function should usually be 
smaller. 


Large negative reflectance errors (case 1) indicate that the algorithm will return 
negative reflectances, if the actual reflectance is weak. The large negative percentage 
reflectance error corresponds to an error magnitude of Ap = -0.012. If the simulation 
surface reflectance was less than +0.012, the algorithm would return a negative reflectance. 

The large errors associated with these uncertainties in the aerosol single scattering 
phase function indicate that the algorithm in its present form is limited to applications where 
the aerosol scattering phase function closely matches the function used in the algorithm. 
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Figure 5. Aerosol single scattering phase functions derived from: 1) bimodal log-normal 

aerosol size distribution used in the unperturbed model ( ), 2) single log-normal 

aerosol size distribution with an index of refraction n = 1.349 - 0.008i, geometric mean 
mass radius r m = 0.6 |im, standard deviation of the logarithm of the radius a = 0.61 (- -), 
and 3) power law aerosol size distribution with v = 4 and an index of refraction 
n = 1.53 - 10* 7 i ( ). 
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However, the algorithm could be easily modified for use with other aerosol size 
distributions by replacing the look-up tables of path radiance, downward flux, 
transmission, and atmospheric backscattering ratio with values generated for a more 
appropriate aerosol size distribution. It is not necessary for these values to be generated by 
the same Dave (1972 a,b,c,d) radiative transfer code used to derive the values used here as 
long as the values tabulated are defined as in eq. (6-8). 

6.1.2 Polarization 

The algorithm which performs the atmospheric corrections uses tabulated radiances 
and fluxes computed with a Dave radiative transfer code. This code is based upon the scalar 
form of the radiative transfer equation, where the assumption is that the light scattered by 
the atmosphere and the earth's surface is unpolarized. If the atmospheric optical thickness 
is small this assumption is satisfactory, since the primary source of the light scattered by 
the atmosphere is the direct, unpolarized sunlight. The polarization errors are caused by 
second and higher order scattering of light that is polarized, and increase with optical 
thickness and decreasing wavelength. Light scattered by the surface, and especially by the 
atmosphere, is polarized. 

Therefore, sensitivity tests are performed to estimate surface reflectance errors 
caused by neglect of polarization. In these tests, polarization is accounted for by another 
Dave code, similar in all other aspects to the scalar code. These sensitivity tests are 
performed for three wavelengths (cases 3,4,5): 0.486 |im (TM band 1), 0.639 pm 
(AVHRR band 1), and 0.845 pm (AVHRR band 2). Errors in derived surface reflectance 
associated with the assumption of unpolarized light are less than 3% for X ^ 0.639 pm 
(cases 4 and 5), but increase to approximately 10% for X = 0.486 pm (TM band 1). 
Although the results are not shown in the tables, the errors associated with the assumption 
of unpolarized light generally decrease with decreasing aerosol optical thickness because 
less multiple scattering occurs. 

6.1.3 Water Vapor and Ozone Absorption 

Another potential source of error which is examined is the effect of using an incorrect 
value of the gaseous absorption. Since most of the gaseous absorption in the Landsat TM 
and NOAA AVHRR visible and near IR bands is contributed by either water vapor or 
ozone, the sensitivity tests are run at those wavelengths where absorption by these 


68 



constituents is the largest. For water vapor this occurs for AVHRR band 2 while for ozone 
this occurs for TM band 2. 

Since the amount of an absorbing gas is variable, the algorithm uses a weighted 
average of the gaseous absorption values computed for the tropical, mid-latitude summer, 
and mid-latitude winter models given in the LOWTRAN code. Most of the weight is given 
to the mid-latitude summer profile (see eq. 16). In the case of AVHRR band 2, the 
algorithm uses a water vapor gaseous absorption optical thickness of Xg = 0.0933. The 
gaseous absorption optical thicknesses computed using the mid-latitude winter and tropical 
profiles are used for the sensitivity studies. Therefore, two tests are made to estimate the 
uncertainty due to water vapor absorption; the first test uses the mid-latitude winter value of 
tg = 0.0486 (case 6) while the second test uses the tropical value of Xg = 0.1235 (case 
7). 


A change in the total amount of water vapor between the algorithm and mid-latitude 
winter profiles (case 6) causes errors as large as 1 1% in the derived surface reflectance for 
AVHRR band 2. Since the difference in the total water vapor amounts between the 
algorithm and tropical profiles is smaller than the difference between the algorithm and mid- 
latitude winter profiles, the magnitude of the error that results from using the tropical 
profile is less— 7% (case 7). These errors can be reduced significantly by applying the water 
vapor correction given in Section 3.3. 

A second set of tests is conducted to determine the sensitivity of the model to changes 
in ozone absorption. Since TM band 2 has the largest ozone absorption, the ozone 
absorption sensitivity tests are run at this wavelength. As in the case of water vapor 
absorption discussed above, the algorithm uses a weighted average of the mid-latitude 
winter, mid-latitude summer, and tropical profiles computed using the LOWTRAN code to 
derive a value of ozone absorption optical thickness of Xg = 0.032 for TM band 2. The 
first of two sensitivity tests are run using the mid-latitude winter ozone absorption optical 
thickness Xg = 0.039 (case 8) while the second test is run using the tropical value of 
Xg = 0.024 (case 9). The magnitude of errors in the the derived surface reflectance in TM 
band 2 caused by variations in mid-latitude winter and tropical profiles to the algorithm 
ozone amounts is 2-4% (cases 8 and 9). 
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6.1.4 Aeroso l Absorption 


The sensitivity of the correction algorithm to aerosol absorption is tested by changing 
the aerosol single scattering albedo ©o, which is the ratio of aerosol scattering to extinction. 
The algorithm uses values of too derived for the 70%-relative-humidity, rural model of 
Shettle and Fenn (1979); in this model the aerosol single scattering albedo varies from 0.95 
atX = 0.486 Jim to 0.86 at X = 2.550 pm. Values of the single scattering albedo for the 
visible spectrum, reported by Waggoner et al. (1981), range from 0.54 <, ©o ^ 0.61 for 
urban industrial areas, 0.73 <, coo ^ 0.87 for urban residential areas, and 
0.89 ^ ©o ^ 1-00 for rural areas. These values agree with the values used by Shettle and 
Fenn. In the sensitivity test (case 10), a value of ©o = 0.84 at 0.486 pm is based upon the 
average measured values in residential urban areas of Michigan and Missouri (Waggoner 
et. al., 1981). Since the algorithm is designed to correct radiances measured over rural 
sites, this should represent a rather extreme (but possible) departure from the usual case. 

The derived surface reflectance is somewhat sensitive to changes in aerosol 
absorption (case 10). This sensitivity appears to depend strongly on geometry, as the 
resulting errors in the derived reflectance vary between -10 to 0% for the case when the 
single scattering albedo decreased from 0.95 to 0.84. 

Another study of the effect of aerosol absorption error on the derived surface 
reflectance is given out of sequence (case 10a ), because the unperturbed model is different. 
The emphasis is on a large solar zenith angle (0 O = 64°), long path through the atmosphere 
(0 = 56°), and rather large aerosol optical thickness (t a = 0.35). Interpolations are made 
on all these variables plus the azimuth. Also, the simulated surface reflectance is high 
(0.6). The error (-0.1 1%) is significant. The error would be less with smaller optical 
thickness, however. 

6.1.5 Vertica l Distribution of Aerosols 

The next sensitivity test deals with errors which result from uncertainty in the vertical 
distribution of aerosols. The altitude distribution of aerosols used in the algorithm is based 
on the 'average' distribution described by Braslau and Dave (1973) which is shown in 
Figure 4. The sensitivity of the results to this distribution is tested by using three different 
profiles which are constructed using the aerosol models described by Shettle and Fenn 
(1979). In these models, the atmosphere is divided into four regions: boundary layer (0 - 
2 km), upper troposphere (2 - 9 km), stratosphere (9 - 30 km), and upper atmosphere 
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(30 - 100 km). All three of the models use the background aerosol models for the 
stratosphere and upper atmosphere regions; the perturbation models differ in the boundary 
layer and upper troposphere regions. The first model (case 11) uses the 50-km- visibility, 
rural model in the boundary layer and the 50-km-visibility, spring/summer model in the 
upper troposphere. The second model (case 12) uses the 23-km-visibility, rural model in 
the boundary layer and the 23-km-visibility, spring/summer model in the upper tropo- 
sphere. The third model (case 13) uses the 10-km- visibility, rural model in the boundary 
layer and the 23-km-visibility, spring/summer model in the upper troposphere. The vertical 
distributions of aerosol corresponding to these three models are shown in Figure 4. These 
profiles are normalized to produce an optical thickness of x a = 0.25 at 0.550 pm. Since 
the average visibility in the midwestem U.S. is about 20 km (Husar and Holloway, 1984), 
a maximum visibility of 50 km and a minimum visibility of 10 km are chosen to represent 
the extremes used for the construction of the altitude distribution of aerosols. 

Uncertainty in the height distribution of aerosols has a negligible effect on the surface 
reflectance when measured from the top of the atmosphere, since the aerosol distribution is 
normalized such that the total aerosol optical thickness below the sensor remains the same 
(cases 1 1, 12 and 13). However, if a sensor is within the atmosphere, there will be small 
differences in the optical thickness of aerosols above and below the sensor, depending 
upon which aerosol distribution is used. The magnitude of the derived reflectance errors is 
less than 7% at aircraft heights (Table 10). 

6.1.6 Bidirectional Reflectance Errors 

The final model error concerns the assumption that the surface is Lambertian. Surfaces 
do not reflect light according to Lambert’s Law as is assumed for the current atmospheric 
correction algorithm. Lee and Kaufman (1986) calculated the error in the derived surface 
reflectance for a model which also assumed Lambert reflection. Their study utilized actual 
surface bidirectional reflectances for pasture as measured by Kriebel (1977). The absolute 
errors in estimates of surface reflectance are a few hundredths when the solar zenith angle 
is small because the surface is nearly Lambertian. The errors are also small where the 
surface reflectance is weak. The derived surface reflectance errors become large (about 
0.1), however, for moderate haze, when the surface reflectance is both high and strongly 
anisotropic. Errors in the derived surface reflectance and radiance, but not the irradiance, 
will not be significant if the surface reflectance is weak; however, the errors can be 
important for strong bidirectional reflectance. 
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6.2 Interpolation Errors 


The correction algorithm is based on a series of look-up tables relating surface 
reflectance to measured upward radiance for various aerosol optical thicknesses and 
geometries. As a result, the algorithm must interpolate to determine the reflectance 
corresponding to an arbitrary input radiance, aerosol optical thickness, gaseous absorption 
optical thickness, altitude, and geometry. (No extrapolation is permitted.) The simulations 
are chosen to show the largest interpolation errors. The interpolation errors are computed 
by comparing the reflectance errors derived from simulated radiances with those derived 
with the correction algorithm. 

The first interpolation sensitivity test (case 14) is performed using the standard test 
input parameters described earlier (X = 0.486 pm, 0© = 30°, 0 = 0°, <|> = 0°), except 
with an aerosol optical thickness of x a = 0.375, which requires the algorithm to interpolate 
between T a = 0.25 and x a = 0.50. Note that the reflectance errors are deviations from the 
correct value of p = 0.05 for cases 14-21 and 24-32; the unperturbed reflectance changes 
only for cases 22 and 23. The next run (case 15) tests the interpolation on solar zenith 
angle. The solar zenith angle interpolation test is run using the standard test input 
parameters except with a solar zenith angle of 26° which requires the algorithm to 
interpolate between 0 O = 20° and 0 O = 30°. Similar tests are run for the observation 
zenith angle (case 16) and the azimuth angle (case 17). The observation zenith angle 
interpolation test was run using 0 = 20°* which requires the algorithm to interpolate 
between 0=18° and 0 = 24°. The azimuth angle interpolation test was run using 
<(> = 145° so that the algorithm interpolated between <|> = 140° and <|> = 150°. Separate 
interpolations on x a , 9, and <j> result in errors less than 1% in the derived surface 
reflectance, while the interpolation on 0 O can result in an error as large as 4 %. 

Errors associated with interpolation on wavelength are shown in cases 18-20. The 
input wavelengths are 0.5, 0.6 and 0.775 pm. In these cases, the error in derived surface 
reflectance is at most 6 %. 

Since the algorithm may also interpolate on gaseous absorption, sensitivity tests are 
also made to determine the uncertainty involved with this interpolation. In cases 21-23 the 
algorithm uses the input gaseous absorption value x g computed using the LOWTRAN 6 
code and the radiances associated with TM band 4 (0.840 pm) (which has a narrow 
bandwidth and relatively small gaseous absorption) to determine the surface reflectance 
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seen by the AVHRR band 2 (which has a wide bandwidth and relatively large gaseous 
absorption). These tests show the algorithm underestimates the surface reflectance by only 
about 1-2%. 

Sensitivity tests are also made to determine the errors associated with adjustments to 
the look-up tables corresponding to different input surface heights. As discussed in section 
3.5, the algorithm makes the adjustment by changing the wavelength slightly (eq. 20). 
Cases 24 and 25 show that this adjustment results in rather large errors in the derived 
surface reflectance of -8%. 

Errors associated with interpolation on the measurement altitude Z are studied. The 
radiative transfer computations are tabulated at three altitudes: 0.45 km, 4.5 km, and 80 
km. In cases 26 to 28 input observation altitudes of Z = 0.23 km, 2.96 km, and 6.5 km are 
used. The derived surface reflectance errors are as large as 8%. 

Finally, errors associated with interpolation for large solar and viewing zenith angles 
are given. Interpolations are made with respect to x a , 6o> 6, and the surface reflectance 
p = 0.05 (cases 29 -32). Case 29 shows that the error is appreciable for a rather large 
optical thickness. This large error depends on the nonlinear change in radiance as the 
azimuth increases from 10° to 20°. If the same test is made, except that the azimuth is 
165°, the error reduces to 8% (case 30). When the optical thickness decreases to x a = 0.1 
(case 31), however, the error has a large negative value. The derived value is p = 0.034, 
compared with the simulation value of p = 0.050. Again, if the azimuth is changed from 
near-forward (case 31) to near- backward (case 32), the magnitude of the error decreases to 
6 %. 

6.3 Input Errors 

A third set of sensitivity tests is performed to estimate the errors in the derived 
surface reflectance resulting from errors or uncertainties in the input data. Each test is 
performed with the algorithm model, except that an incorrect value of one of the input 
parameters is used. 

The first test introduces a 5% error into the input radiance (case 33). As would be 
expected, the algorithm is very sensitive to such an error. If the input radiance error 
increases to 10% (case 34 in Tables 8 and 10), the reflectance error doubles to 28%. The 
relative errors in the derived surface reflectance can be as large as two to three times the 
relative errors in the input radiance. 
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Errors in the input optical thickness occur because of measurement errors, and the 
optical thickness is not measured at the same time and place of the remote measurement. 
Tests are made with an error in the input aerosol optical thickness of Ata = 0.05 for the 
case where the correct input is t a = 0.25 (case 36), and for At a = 0.1 when T a = 1.0 (case 
38). The derived surface reflectance errors caused by optical thickness errors are largest 
when the surface reflectance is weak. Hence, the errors are computed for visible reflectance 
of visible light from vegetation (p = 0.05). The resulting surface reflectance error of 
0.003 is insensitive to errors of 0.05 for moderate aerosol optical thickness (case 36). For 
t a = 1.0 and larger aerosol optical thickness error of 0.1, the error in the derived surface 
reflectance increases to Ap = 0.014 (case 38). 

The solar zenith angle test is run by introducing an error of A0o = 2° (case 40) when 
the correct input is 0 O = 30°. Similarly, the sensitivity of the results to errors in the input 
observation zenith angle is tested by introducing an error of A0 = 2° (case 41) when the 
correct input is 0 = 0°. Finally, the sensitivity to errors in the input azimuth angle is tested 
by introducing an error of A<|> = 2° for the case when the correct input is <J> = 0° (case 
44). Errors in the input values of solar zenith, observation zenith, and azimuth angles 
generally result in negligible errors of less than 2% in the derived surface reflectance. 

7.0 Conclusion 

An algorithm is developed to account for atmospheric effects when deriving surface 
reflectance properties from visible and near-infrared radiances measured by aircraft or 
satellite over rural areas. The radiance that would be measured for a given surface 
reflectance can be derived, also. The algorithm uses a tabulated set of radiances computed 
for various wavelengths, solar and observation angles, and aerosol optical thicknesses. All 
aerosol parameters have been assumed, except for the aerosol optical thickness, which is an 
input value. Since the algorithm performs essentially interpolations, it is fast; therefore, it 
is well suited for reducing observations in many wavelengths. Otherwise, the effect of the 
atmosphere requires many radiative transfer computations. 

Large errors in derived parameters, rather than rms errors are estimated. Among the 
largest model errors are those caused by uncertainties in the aerosol scattering phase 
function; in this case surface reflectance and radiance errors reach ±20-30%. Thus, in its 
current configuration, the algorithm is suitable for only a rural, bimodal log-normal aerosol 
size distribution. However, the algorithm could be easily modified for use with other 
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aerosol size distributions by replacing the radiance look-up tables with values generated for 
a more appropriate aerosol size distribution. 

Errors in the derived surface reflectance can be large when either the slant path through 
the atmosphere of sunlight or light reflected from the ground is long.The uncertainty in the 
amount of water vapor causes an error of 5% in the reflectance for AVHRR band 2. This 
error can be reduced significantly by using measurements of the total amount of water 
vapor at the time of measurement The mesh for the tabulated radiation parameters is fine 
enough so that linear interpolation results in small reflectance errors (<4%). 

Of the errors in the input parameters required by the correction algorithm, mors in the 
measured radiance can result in large errors in the derived radiance and reflectance but not 
irradiance. Absolute measured radiance errors of 5-10% are expected. Therefore, the 
derived surface radiances and reflectances will deviate by at least 5-10% from the correct 
values even without any errors in the atmospheric correction algorithm. Large errors in the 
derived surface reflectance can also result from errors in the input values of aerosol optical 
thickness. Optical thickness errors should be less than 0.05 so that the corresponding 
surface reflectance errors are less than 10% for a dark surface ( p = 0.05). 

Acknowledgement: We appreciate Mr. Brian Markham’s early efforts to use the tables, and 
thereby bring to our attention places to improve their accuracy. 
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9.0 FORTRAN Listing 

The FORTRAN code for the atmospheric correction algorithm is given in the pages 

which follow. The line numbers are listed to the right of each line. 
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THIS SOFTWARE WAS DEVELOPED BY FIF00010 

SHANA MATTOO, ARC AND GSFC FIF00020 

C109D, BLD22, TEL# 286-2120 FIF00030 

********************************************************************* F IF 00040 

* THIS MAIN PROGRAM IS BASED ON THE ATMOSPHERIC CORRECTION OF UPWARD* FIF00050 

* RADIANCE MEASURED FROM AIRCRAFT OR SATELLITE . IT DERIVES DOWNWARD * FIF00060 

* IRRADIANCE, UPWARD RADIANCE, AND REFLECTANCE ALL FOR THE GROUND * FIF00070 

* LEVEL. THIS PROGRAM ALSO HAS A OPTION TO COMPUTE THE UPWARD * FIF00080 

* RADIANCE IF THE SURFACE REFLECTANCE IS PROVIDED. * FIF00090 

********************************************************************* FIF00100 


SUBROUTINES USED: 

1. READ IN 

2. FINDW 

3 . INTSFX 

4. INTHGH 

5 . INTERP 

6 . INTEXP 

7. SYSTEM SUBROUTINE IMSL(ZXGSN) (LOOK IN SUBROUTINE INTEXP) 
(BE SURE TO LINKTO PROPER SYSTEMS LIBRARY CONTAING IMSL) 


FIF00110 

FIF00120 

FIF00130 

FIF00140 

FIF00150 

FIF00160 

FIF00170 

FIF00180 

FIF00190 

FIF00200 

FIF00210 

FIF00220 

FIF00230 

FIF00240 

FIF00250 

FIF00260 

FIF00270 

FIF00280 

FIF00290 


REAL * 4 MWAV, MTAU, MINT, MTHET0 , MTHET, MPHI , MHGHT, RHOS, MRHO, SBARN 
REAL *4 HGHT (3) , THE (13) , INT (9, 4, 13, 19) , WAV, OPTH (4) , FDOWN (9, 4) 

REAL *4 SBAR (4),PIT(13,4), THET0 ( 9) , NEW (30) , YY ( 1 ) , NEWN (30 ) 

REAL *4 SBARW (2,4), PITW (2,13,4,3), FDOWNW (2 , 9, 4 ) 

REAL *4 INTW (2,9,4,3,13,19), INTWH (2, 9, 4, 13, 19) , PITWH (2,13,4) 

REAL *4 NEWINT (4, 13, 19) ,NFDOWN(4) ,NEWNEW(4, 13) ,FINT(4) , FT (4) 

REAL *4 F, AIRR, ARAD, PHI (19) 

REAL *4 F0IRR(100) , WAVE (100) ,R(365) 

REAL *4 LEAP ( 12 )/30., 29. ,31. ,30. ,31. ,30. ,31. ,31. ,30. ,31. ,30. ,31./ FIF00300 
REAL *4 NLEAP (12) /30 . , 28 . , 31 . , 30 . , 31 . , 30 . , 31 . , 31 . , 30 . , 31 . , 30 . , 31 . /FIF00310 
REAL *4 WAVEN(8) / .486, .587, .639, .663, .837, .845,1.663,2.189/ FIF00320 

REAL *4 TAUGH(8) /. 0070, .0320, .0247, .0174, .0021, .0152, .0077, .0091/ FIF00330 
REAL *4 TAUGL(8)/. 0130, .0142, .0208, .0218, .0628, .1156, .1380, .0930/ FIF00340 
CHARACTER * 1 LINE (132) 

CHARACTER * 1 LINE2(80) 

CHARACTER * 1 LL(30) 

INTEGER ANGLE (19) 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


DEFINE PI 

PI=ARCOS (-1.00) 

********************************************* ************** ********* 
* CALL SUBROUTINE READIN * 

******************************************************************** 
SUBROUTINE READIN READS THE INPUT DATA AND WRITES IT ON UNIT 6, TO 
PERFORM A CHECK ON INPUT DATA. 

CALL READIN 


LINE2 READS THE TITLE FOR THE VALUES OF FOLLOWING INPUT QUANTITIES . 

READ THE IOPT, MOPT AND NOPT. IOPT IS OPTION FOR UNITS OF RADIANCE. 

IF OPTION IS 1, UNITS SHOULD BE IN REFLECTANCE UNITS; IF IOPT=2 
THE UNITS OF REFLECTANCE SHOULD BE IN ABSOLUTE UNITS . IF MOPT IS 1 
THE DAY, MONTH AND YEAR FOR WHICH DATA IS MEASURED SHOULD BE 
GIVEN; AND IF MOPT IS 2 THE DAY OF THE YEAR AND THE YEAR SHOULD BE 
GIVEN .NOPT IS OPTION FOR COMPUTING SURFACE REFLECTANCE OR RADIANCE 
IF NOPT IS 1 SURFACE REFLECTANCE IS COMPUTED AND IF NOPT IS 2 RADIANCEFIF00580 
IS COMPUTED. FIF00590 

FIF00600 
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READ(5,1003)LINE2 
READ (5, 9) IOPT, MOPT, NOPT 


C 

C LINE2 READS THE TITLE FOR THE VALUES OF FOLLOWING INPUT QUANTITIES: 

C INUM THE NUMBER OF DATA POINTS TO BE PROCESSED, MONTH, DAY, YEAR AND 
C TIME ZONE FOR WHICH MEASURED DATA ARE PROVIDED. 

C 

READ (5, 1003) LINE2 

IF (MOPT .EQ. 1) READ (5, 10) INUM, IMONTH, IDAY, I YEAR, LL 
IF (MOPT .EQ. 2) READ (5, 10) INUM, IDAY, IYEAR, LL 
READ (5,1003) LINE2 
C 

C WRITE THE LABELS FOR THE OUTPUT. 

C 

WRITE (56, 1108) 

WRITE (56, 1109) 

IF (MOPT . EQ . 1 ) WRITE (6,11) IMONTH, IDAY, IYEAR, LL 
IF (MOPT .EQ . 1) WRITE (56, 11) IMONTH, IDAY, IYEAR, LL 
IF (MOPT . EQ . 2 ) WRITE (6,21) IDAY, IYEAR, LL 
IF (MOPT .EQ. 2) WRITE (56,21) IDAY, IYEAR, LL 
WRITE (6, 12) 

WRITE (56,22) 

C 

CREAD LABS AND NECKEL DATA FOR SOLAR FLUX 
C 

DO 121 IJ - 1,4 

121 READ (20, 1003) LINE 
IK*=1 

DO 122 IJ - 1,20 

READ (20,23) WAVE ( IK) ,F0IRR(IK) ,WAVE(IK+1) ,F0IRR(IK+1) ,WAVE(IK+2) , 

1 F0IRR (IK+2) 

IK=IK+3 

122 CONTINUE 
IK=IK-1 

C 

C START OF LOOP FOR THE NUMBER OF THE DATA POINTS TO BE PROCESSED 
C LFILE IS USED TO SET THE VALUE FOR THE FILE NUMBER NFILE TO READ 
C THE LOOK-UP TABLE FOR CHOSEN WAVELENGTH. 

C 

LFILE=6 

2000 DO 999 NUM“1, INUM 

IF (NUM .GT.l) REWIND NFILE1 
IF (NUM .GT.l) REWIND NFILE2 

q ******************************************************************** 

C * READ MEASURED PARAMETERS * 

q ******************************************************************** 

C TIME IN CDT (MTIME IN HOURS AND MINUTES) , WAVELENGTH (MWAV IN 
C MICROMETERS) , OPTICAL THICKNESS (MTAU) FOR MWAV, SOLAR ZENITH ANGLE 
C MTHET0 IN DEGREES), OBSERVATION SCAN ANGLE (MTHET IN DEGREES), 

C OBSERVATION AZIMUTH ANGLE (MPHI IN DEGREES), OBSERVATION HEIGHT (MHGHT 
C IN KILOMETERS) , INTENSITY (MINT IN REFLECTANCE OR ABSOLUTE UNITS) IF 
C NOPT IS 1; AND SURFACE REFLECTANCE (MRHO) IF NOPT IS 2. 

C TAUGHW (ABSORPTION FOR CARBONDIOXIDE AND OZONE) , TAUGLW (ABSORPTION 
C FOR WATER) FOR MEASURED WAVELENGTH, HEIGHT FROM THE SURFACE ABOVE SEA- 
C LEVEL (SHGHT IN KILOMETERS) . (LOOK DOCUMENTATION PAGE 19) 

C 

IF (NOPT . EQ . 1 ) 

1READ (5,13) MTIME, MWAV, MTAU, MTHET 0, MTHET, MPHI, MHGHT, MINT, TAUGLW, 

1 TAUGHW, SHGHT 
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IF (N0PT.EQ.2) 

1READ (5, 13) MT IME , MWAV, MTAU , MTHETO , MTHET , MPH I , MHGHT , MRHO , TAUGLW, 

1 TAUGHW, SHGHT 

C 

C IF TAUGLW (ABSORPTION FOR WATER) IS 0.0 THE DEFAULT VALUE IS COMPUTED 
C BY INTERPOLATING LINEARLY FOR THE MEASURED WAVELENGTH MWAV. IF MWAV IS 
C OUT OF RANGE NO EXTRAPOLATION IS ALLOWED, THE PROGRAM SENDS CONTROL TO 
C PROCESS NEW DATA POINT. 

C 

IF (TAUGLW .GT. 0.0) THEN 

GO TO 360 

ELSE 

CALL INTERP (1,8, MWAV, WAVEN, TAUGL, YY, I JK, NUM, LOPT) 

IF (LOPT . EQ . 0 ) THEN 
WRITE ( 6 , 2 4 ) MWAV, NUM 
GO TO 999 
ELSE 

TAUGLW=YY ( 1 ) 

ENDIF 

END IF 
C 

C IF TAUGHW (ABSORPTION FOR GASES) IS 0.0 THE DEFAULT VALUE IS COMPUTED 
C BY INTERPOLATING LINEARLY FOR THE MEASURED WAVELENGTH MWAV. IF MWAV IS 
C OUT OF RANGE NO EXTRAPOLATION IS ALLOWED, THE PROGRAM SENDS CONTROL TO 
C PROCESS NEW DATA POINT. 

C 

360 CONTINUE 

IF (TAUGHW .GT. 0.0) THEN 

GO TO 350 

ELSE 

CALL INTERP (1,8, MWAV, WAVEN , TAUGH , YY , I JK , NUM , LOP T ) 

IF (LOPT . EQ . 0 ) THEN 
WRITE ( 6 , 2 5 ) MWAV, NUM 
GO TO 999 
ELSE 

TAUGHW=YY ( 1 ) 

ENDIF 

ENDIF 

350 CONTINUE 
C 

q** ************************************************* ****************** 
C* CALL TO SUBROUTINE FINDW * 
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Q* *********************************************** ********************* 

C SUBROUTINE FINDW PICKS TWO WAVELENGTHS WAV1 AND WAV2 BETWEEN WHICH 
C MEASURED WAVELENGTH MWAV LIES AND SETS FILE NUMBERS NFILE1 AND NFILE2 
C TO READ THE REQUIRED DATA SET FOR THE TWO WAVELENGTHS FROM LOOK-UP 
C TABLE. THE MWAV SHOULD BE IN RANGE OF 0.486 - 2.189 UM, NO EXTRAPOLAT- 
C ION IS ALLOWED. 

C 

CALL FINDW (MWAV, I I , NUM, KOPT , IMM, WAV1 , WAV2 ) 

IF(KOPT.EQ.O) GO TO 999 
NFILE1 = LFILE+II 
NFILE2 = LFILE+IMM 

C 

C ******************************************************************** 

C * CALL TO SUBROUTINE INTSFX * 

C ******************************************************************** 

C SUBROUTINE INTSFX INTERPOLATES TO COMPUTE THE VALUES OF R FOR 
C 365 DAYS. THESE VALUES WILL BE USED TO COMPUTE THE CORRECTED VALUES 
C OF SOLAR FLUX FOR THE DAY MONTH AND YEAR FOR WHICH MEASURED DATA 
C IS GIVEN. R IS SUN-EARTH DISTANCE . 
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c 

CALL INTSFX(R) 

CHANGE MONTH AND DAY TO THE JULIAN DATE IF MOPT IS 1. 

IF (MOPT. EQ. 2) THEN 
IQ=IDAY 
ELSE 

IS=IDAY+1 
IM= IMONT H - 1 
IQ=0 

DO 20 IJ=1,IM 
IF (MOD ( I YEAR, 4 ) . EQ . 0 ) THEN 
IQ-IQ+LEAP (IJ) 

ELSE 

IQ=IQ+NLEAP (IJ) 

ENDIF 

20 CONTINUE 
IQ=IQ+IS 
ENDIF 

COMPUTE SOLAR FLUX FOR THE MEASURED WAVELENGTH MWAV FROM LABS AND 
NECKEL DATA BY INTERPOLATION 

IJK-0 
LNUM-0 

CALL INTERP <1, IK, MWAV, WAVE, F0 IRR, YY, I JK, LNUM, LOPT) 


COMPUTE CORRECTED SOLAR FLUX FOR THE DAY OF YEAR FOR WHICH 
MEASUREMENT ARE MADE. (SEE DOCUMENTATION EQ.6) 

R(IQ)-1. 

FFLUX-YY(l) / (R(IQ) *R(IQ) ) 

C 

C 

C COMPUTE THE OBSERVATION ZENITH ANGLE FOR THE SCAN ANGLE OF 
C THE SATELLITE OR AIRCRAFT (SEE DOCUMENTATION EQ.9) 

C RAD=RADIUS OF EARTH 
C 

DEGRAD=ARCOS (-1 . ) /180 . 

RADDEG=180 . /ARCOS (-1 . ) 

RAD=6370 . 

MTHET=ASIN ( ( 1 . + (MHGHT/RAD) ) * (SIN (MTHET*DEGRAD) ) ) *RADDEG 
C 

C IF IOPT EQ 1 CONVERT MEASURED RADIANCE TO ABSOULTE UNITS. 

C IF IOPT IS 2 CONVERT MEASURED RADIANCE FROM ABSOLUTE UNITS TO 
C REFLECTANCE UNITS. 

C 

AMUOCOS (MTHET0*DEGRAD) 

IF ( IOPT . EQ . 1 .AND. NOPT.EQ.l ) AMINT- (MINT*FFLUX*AMUO ) /PI 
IF (IOPT . EQ. 2 .AND. NOPT.EQ.l) AMINT=MINT 

IF ( IOPT . EQ.2 .AND. NOPT.EQ.l) MINT= (MINT*PI) / (FFLUX*AMU0) 
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Q ******************************************************************** 

C * READ FROM LOOK-UP TABLE. * 

Q ******************************************************************** 

Q ******************************************************************** 

C * WAV=WAVELENGTHS (8) (.486, .587, .639, .663, .837, .845, 1.663,2.189) * 

C * MICROMETERS. * 

C * LTHETO “NUMBER OF THETAO (SOLAR ZENITH ANGLE) * 

C * THET0= (10, 20,30, 40, 50, 60, 66,72, 78 DEGREES) * 

C * LTHE “NUMBER OF THE (OBSERVATION ZENITH ANGLE) * 

C * THE= (0,6, 72 DEGREES EVERY 6 DEGREES) * 

C * LPHI “NUMBER OF PHI (OBSERVATION AZIMUTH ANGLE) * 

C * PHI=(0 r 5, 10, 20,30, 40, 60, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, * 

C * 175,180 DEGREES) * 

C * LTAU “NUMBER OF OPTH (OPTICAL THICKNESS) * 

C * OPTH=(0.0, .25, .50,1.0) * 

C * LHGHT - NUMBER OF HGHT (OBSERVATION HEIGHTS) * 

C * HGHT= (80.0,4.5, .45 KM) * 

0 ******************************************************************** 
LTHET0=9 
LTHE=13 
LPHI“19 
LTAU=4 
LHGHT“3 

C READ OBSERVATION ZENITH ANGLE (THE) , AND OBSERVATION AZIMUTH ANGLE 
C FROM LOOK-UP TABLE. STORE INTEGER ARRAY ANGLE IN REAL ARRAY PHI . 

C 

READ (NFILE1, 1000) (THE (I) , I-1,LTHE) 

READ (NFILE1, 1001) (ANGLE (I) , I“1,LPHI) 

READ (NFILE2, 1000) (THE (I) , 1=1, LTHE) 

READ (NFILE2, 1001) (ANGLE (I) , 1-1, LPHI) 

DO 30 IPHI=1 , LPHI 
30 PHI (IPHI) -ANGLE (IPHI) 

C 

C READ WAVELENGTH (WAV) , OPTICAL THICKNESS (OPTH) , SOLAR ZENITH ANGLE 
C (THET0), REFLECTANCE OF ATMOSPHERE (SBARW) , FLUX INCIDENT ON 
C SURFACE (FDOWNW) FROM LOOK-UP TABLE. 

C 

C IF WAVELENGTH IS 1.663 OR 2. 189, THE LOOK-UP TABLE IS AVAILABLE ONLY 
C FOR 2 OPTICAL THICKNESSES OF 0 . 0 AND .25. 

C 

DO 9995 IWAV =1,2 
IF ( IWAV .EQ. 1) NFILE = NFILE1 
IF (IWAV .EQ. 2) NFILE - NFILE2 
IF (NFILE .GT. 12)LTAU=2 
DO 100 ITHET0=1, LTHET0 
DO 101 ITAU=1, LTAU 
DO 102 IHGHT=1, LHGHT 

READ (NFILE, 1002 ) WAV, OPTH ( ITAU) , THET0 (ITHET0) ,SBARW(IWAV, ITAU) , 

1 FDOWNW ( IWAV, ITHET0 , ITAU) 

READ (NFILE, 1003) LINE 
C 

C READ ATMOSPHERIC RADIANCE (INTW) AS FUNCTION OF SOLAR ZENITH ANGLE, 

C OPTICAL THICKNESS, HEIGHT, OBSERVATION ZENITH ANGLE, AND OBSERVATION 
C AZIMUTH ANGLE FROM LOOK-UP TABLE. 

C 

DO 103 ITHE=1, LTHE 

READ (NFILE, 1004) (INTW (IWAV, ITHET0, ITAU, I HGHT, ITHE, IPHI) , IPHI=1, 10) 
READ (NFILE, 1004) (INTW(IWAV, ITHET0, ITAU, I HGHT, ITHE, IPHI) , 

1IPHI=11, 19) 

103 CONTINUE 
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c 

C READ TRANSMISSION FROM SURFACE (PITW) AS FUNCTION OF OBSERVATION 
C ZENITH ANGLE, OPTICAL THICKNESS, AND HEIGHT. 

C 

READ (NFILE, 1005) (PITW(IWAV, ITHE, ITAU, IHGHT) , ITHE=1, LTHE) 

102 CONTINUE 
101 CONTINUE 
100 CONTINUE 
9995 CONTINUE 
C 

0 ******************************************************************** 
C *CALL TO SUBROUTINE INTHGH * 

Q ******************************************************************** 

C SUBROUTINE INTHGH INTERPOLATES THE VALUES OF INTW, AND PITW FOR THE 
C MEASURED HEIGHT MHGHT (SEE DOCUMENTATION SECTION 3.6) 

C 

CALL INTHGH (LTAU, LTHE, LPHI, LTHET0 , INTW, MHGHT, PITW, INTWH, PITWH) 

C 

C ADJUST LOOK-UP TABLE FOR FOR SURFACE LEVEL (SEE DOCUMENTATION EQ. 

C 18, 19 AND 20) 

C 

WAV1 - WAV1 * EXP ( ( SHGHT- . 4 ) / 3 6 . ) 

WAV2 - WAV2 * EXP ( ( SHGHT- . 4 ) / 3 6 . ) 

C 

C INTERPOLATE INTWH, PITWH, SBARW, FDOWNW FOR THE MEASURED WAVELENGTH 
C MWAV. (SEE DOCUMENTATION 3.1, PAGE 12 - 15) 

C 

ALPHA - ( ALOG (MWAV) - ALOG(WAV2)) / (ALOG(WAVl) - ALOG(WAV2)) 

DO 200 ITAU - 1 , LTAU 

SBAR(ITAU)- EXP ((ALPHA * ALOG (SBARW ( 1, ITAU) ) ) + 

1 ((1. -ALPHA) * ALOG (SBARW (2, ITAU) )) ) 

DO 201 ITHE - 1, LTHE 

PIT (ITHE, ITAU) - EXP ((ALPHA * ALOG (PITWH (1, ITHE, ITAU) )) + 

1 ((1. -ALPHA) * ALOG (PITWH (2, ITHE, ITAU) )) ) 

DO 202 IPHI - 1, LPHI 
DO 203 ITHET0 - 1, LTHET0 

FDOWN(ITHET0, ITAU) - EXP ( (ALPHA * ALOG (FDOWNW (1, ITHET0, ITAU) ) ) + 
1 ((1. -ALPHA) * ALOG (FDOWNW (2, ITHET0, ITAU) ) ) ) 

INT<ITHET0, ITAU, ITHE, IPHI) - 
1 EXP ((ALPHA * ALOG (INTWH (1,ITHET0, ITAU, ITHE, IPHI) ) ) 

1 + ((1. -ALPHA) * ALOG (INTWH (2, ITHET0, ITAU, ITHE, IPHI) )) ) 

203 CONTINUE 
202 CONTINUE 
201 CONTINUE 
200 CONTINUE 
C 

C COMPUTE DELTAH AND DELTAL AS EXCESS OR DEFICIT IN UPPER (OZONE, CARBO 
C N-DIOXIDE) AND LOWER (WATER) ATMOSPHERE RESPECTIVELY. 

C 

AVTAGH = EXP ((ALPHA * ALOG (TAUGH (II) ) ) + 

1 (d. -ALPHA) * ALOG (TAUGH (IMM) )) ) 

DELTAH = TAUGHW - AVTAGH 
AVTAGL= EXP ((ALPHA * ALOG (TAUGL (II) ) ) + 

1 ((1. -ALPHA) * ALOG (TAUGL (IMM) )) ) 

DELTAL = TAUGLW - AVTAGL 
C 

C COMPUTE FDOWN, SBAR, PIT AND INT AFTER ADJUSTING FOR THE EXCESS AND 
C DEFICIT OF GASEOUS ABSORPTION ( SEE EQ . 17A, 17B, 17C, 17D) 

C 
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DO 908 ITHET0-1,LTHET0 
AMUTHO - COS (THETO (ITHETO) * DEGRAD) 

THDOWN-EXP (-DELTAH/ AMUTHO) 

TLDOWN-EXP (-DELTAL/ AMUTHO) 

DO 908 ITAU-1,LTAU 

908 FDOWN (ITHETO, ITAU)- FDOWN (ITHETO, ITAU) *THDOWN*TLDOWN 
DO 909 ITAU-1,LTAU 

909 SBAR(ITAU)- SBAR(ITAU) * EXP (-DELTAL * 2) 

DO 910 ITAU-1, LTAU 

DO 910 ITHE— 1, LTHE 
AMUTH = COS (THE ( ITHE) * DEGRAD) 

THUP-EXP (-DELTAH/ AMUTH) 

TLUP-EXP (-DELTAL/ AMUTH) 

910 PIT (ITHE, ITAU) -PIT (ITHE, ITAU) *THUP*TLUP 
DO 912 ITAU-1,LTAU 
DO 912 ITHETO— 1, LTHET0 
AMUTHO = COS (THETO (ITHETO) * DEGRAD) 

DO 912 ITHE-1,LTHE 
AMUTH = COS (THE (ITHE) * DEGRAD) 

DO 912 IPHI-1,LPHI 

912 INT (ITHETO, ITAU, ITHE, IPHI) -INT (ITHETO, ITAU, ITHE, IPHI) * 
1 (EXP (- ( (DELTAL/2 . ) +DELTAH) * ( (1 ./AMUTH) + (1 ./AMUTHO) ) ) ) 


Q ★*★**★****★*★★*★*****★*★*★****★★★★***★*★★*****★***★*****★★***★*★★*** 

C * INTERPOLATION OF ATMOSPHERIC RADIANCE AND FLUX INCIDENT ON GROUND* 
C * ON MEASURED SOLAR ZENITH ANGLE. * 

C INT (THETO, OPTH,HGHT, THE, PHI) IS INTERPOLATED FOR 
C MEASURED MTHET0 TO GET NEW VARAIBLE NEWINT (TAU, THE, PHI) . 

C FDOWN (THETO, TAU) IS INTERPOLATED FOR MEASURED MTHET0 TO 
C GET NEW VARIABLE NFDOWN ( TAU ) . 

C 

701 CONTINUE 

DO 240 1-1,30 
NEW (I) -0.0 
240 NEWN ( I ) =0 . 0 

DO 104 ITAU— 1, LTAU 
DO 105 ITHE- 1, LTHE 
DO 106 IPHI— 1, LPHI 
DO 107 ITHETO— 1, LTHET0 

NEW (ITHETO) -INT (ITHETO, ITAU, ITHE, IPHI) 

NEWN (I THETO) -FDOWN (I THETO, ITAU) 

107 CONTINUE 
IJK-1 
C 

C CALL SUBROUTINE INTERP TO INTERPOLATE INT. 

C IF MTHET0 IS OUT OF RANGE OF VALUES OF THETO SUBROUTINE INTERP 
C RETURNS THE VALUE OF LOPT - 0 AND THE NEW DATA POINT WILL BE READ. 

C 

CALL INTERP ( 1 , LTHET0 , MTHET0 , THETO, NEW, YY, I JK, NUM, LOPT) 
IF(LOPT.EQ.O) GO TO 999 
NEWINT ( ITAU, ITHE, IPHI ) -YY ( 1 ) 

IJK-IJK+1 

C 

C CALL SUBROUTINE INTERP TO INTERPOLATE FDOWN. 

C IF MTHET0 IS OUT OF RANGE OF VALUES OF THETO SUBROUTINE INTERP 
C RETURNS THE VALUE OF LOPT = 0 AND THE NEW DATA POINT WILL BE READ. 

C 
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CALL INTERP ( 1 , LTHETO , MTHETO, THETO , NEWN, YY, I JK, NUM, LOPT) 

IF (LOPT.EQ. 0) GO TO 999 
NFDOWN ( ITAU) =YY ( 1 ) 

106 CONTINUE 
105 CONTINUE 
104 CONTINUE 

Q ******************************************************************** 

C * INTERPOLATION OF THE ATMOSPHERIC RADIANCE ON MEASURED AZIMUTH * 
C * ANGLE. * 

c ******************************************************************** 
C NEWINT (TAU, THETA, PHI) IS INTERPOLATED FOR MEASURED MPHI 
C TO GET NEW VARAIABLE NEWNEW (TAU, THETA) . 

C 

DO 241 1=1,30 
NEW (I) =0 . 0 

241 NEWN ( I ) =0 . 0 

DO 108 ITAU=1, LTAU 
DO 109 ITHE=1, LTHE 
DO 110 IPHI=1, LPHI 
NEW (IPHI)=NEWINT( ITAU, ITHE,IPHI) 

110 CONTINUE 
IJK=3 
C 

C CALL SUBROUTINE INTERP TO INTERPOLATE ON PHI. 

C IF MPHI IS OUT OF RANGE OF VALUES OF PHI SUBROUTINE INTERP 
C RETURNS THE VALUE OF LOPT = 0 AND THE NEW DATA POINT WILL BE READ. 

C 

CALL INTERP (1, LPHI, MPHI, PHI, NEW, YY, I JK, NUM, LOPT) 

IF (LOPT.EQ. 0) GO TO 999 
NEWNEW ( ITAU , I THE ) =YY ( 1 ) 

109 CONTINUE 
108 CONTINUE 

C* INTERPOLATION OF ATMOSPHERIC RADIANCE AND TRANSMISSION/PI FROM * 
C* SURFACE ON MEASURED OBSERVATION ZENITH ANGLE. * 

Q**i**********************<r************************t****************** 

C NEWNEW (OPTH, THE) IS INTERPOLATED FOR MEASURED THETA (MTHET) 

C TO GET NEW VARAIABLE FINT (TAU) . TRANSMISSION FACTOR PIT IS 
C INTERPOLATED TO GET FT (TAU) . 

C 

DO 242 1=1,30 
NEW (I) =0.0 

242 NEWN ( I ) =0 . 0 

DO 111 ITAU=1, LTAU 
DO 112 ITHE=1, LTHE 
NEW ( ITHE) =NEWNEW ( ITAU, ITHE) 

NEWN (ITHE) =PIT (ITHE, ITAU) 

112 CONTINUE 
I JK=4 
C 

C CALL SUBROUTINE INTERP TO INTERPOLATE RADIANCE ON MTHET. 

C IF MTHET IS OUT OF RANGE OF VALUES OF THE SUBROUTINE INTERP 
C RETURNS THE VALUE OF LOPT = 0 AND THE NEW DATA POINT WILL BE READ 
C 

CALL INTERP (1, LTHE, MTHET, THE, NEW, YY, IJK, NUM, LOPT) 

IF (LOPT.EQ. 0) GO TO 999 
FINT (ITAU) =YY(1) 

I JK=5 
C 
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C CALL SUBROUTINE INTERP TO INTERPOLATE PIT ON MTHET . 

C IF MTHET IS OUT OF RANGE OF VALUES OF THE SUBROUTINE INTERP 
C RETURNS THE VALUE OF LOPT = 0 AND THE NEW DATA POINT WILL BE READ. 

C 

CALL INTERP ( 1, LTHE, MTHET, THE, NEWN, YY, I JK, NUM, LOPT) 

IF(LOPT.EQ.O) GO TO 999 
FT (ITAU) =YY (1) 

111 CONTINUE 

Q********************************************************************* 

c* * 

C* INTERPOLATE RADIANCE (FINT) , FLUXDOWN (NFDOWN) , TRANSMISSION (FT) AND* 

C* SBAR(SBAR) ON MEASURED OPTICAL THICKNESS (MTAU) . * 

C* * 

q* ******************************************************************** 

C CALL SUBROUTINE INTEXP WHICH PERFORMES EXPONANTIAL INTERPOLATION, 

C IF VALUE OF IGNORE RETURNED FROM SUBROUTINE INTEXP IS 1 LINEAR INTER- 
C POLATION IS PERFORMED USING SUBROUTINE INTERP. 

C 

CALL INTEXP (OPTH, FINT, MTAU, YY, LTAU, IGNORE) 

IF MTAU IS OUT OF RANGE OF VALUES OF OPTH SUBROUTINE INTERP 
RETURNS THE VALUE OF LOPT - 0 AND THE NEW DATA POINT WILL BE READ 


IJK=6 

IF (IGNORE .EQ . 1) THEN 

CALL INTERP (1, LTAU, MTAU, OPTH, FINT, YY, I JK, NUM, LOPT) 

IF ( LOPT . EQ . 0 ) GO TO 999 
END IF 

FINTN=YY(1) 

CALL INTEXP (OPTH, NFDOWN, MTAU, YY, LTAU, IGNORE) 

IF (IGNORE . EQ . 1 ) THEN 

CALL INTERP ( 1 , LTAU, MTAU, OPTH, NFDOWN, YY, I JK, NUM, LOPT) 

IF (LOPT . EQ . 0 ) GO TO 999 
END IF 

FDOWNN=YY ( 1 ) 

CALL INTEXP (OPTH, FT, MTAU, YY, LTAU, IGNORE) 

IF (IGNORE . EQ. 1) THEN 

CALL INTERP ( 1 , LTAU, MTAU, OPTH, FT, YY, I JK, NUM, LOPT) 

IF (LOPT.EQ. 0) GO TO 999 
END IF 
FTN=YY ( 1 ) 

CALL INTEXP (OPTH, SBAR, MTAU, YY, LTAU, IGNORE) 

IF (IGNORE . EQ . 1 ) THEN 

CALL INTERP (1, LTAU, MTAU, OPTH, SBAR, YY, IJK, NUM, LOPT) 

IF (LOPT . EQ . 0 ) GO TO 999 
END IF 

SBARN-YY(l) 

(3********************************************************************* 

C* COMPUTE SURFACE REFLECTANCE RHO (SEE DOCUMENTATION EQ.2,3) * 

C********************************************************************* 

C COMPUTE F = (MINT-FINTN) / ( FDOWNN * FTN ) 

C SURFACE REFLECTANCE RHO=F/ (1+ (F*SBARN (TAU) ) ) 

C IF NOPT IS 2 REFLECTANCE IS COMPUTED. 

C 

IF (NOPT. EQ. 2) GO TO 56 
F= (MINT-FINTN) / (FDOWNN* FTN) 

RHOS=F / ( 1 . + (F*SBARN) ) 
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Q ******************************************************************** 

C * COMPUTE IRRADIANCE IN WATTS /METER* *2 /UM, AND ABSOLUTE RADIANCE * 

C * IN WATTS/METER* *2/UM/SR. * 

Q ******************************************************************** 

C COMPUTE IRRADIANCE AIRR- (FDOWNN*FFLUX* MUO ) / (1-SBARN . RHOS) 

C COMPUTE ABSOLUTE RADIANCE- (RHOS*FDOWN*FFLUX*AMUO) / (PI* (l-SBAR*RHOS) ) 

C 

AIRR= (FDOWNN * FFLUX* AMUO) / ( 1 . -SBARN*RHOS) 

ARAD- (RHOS*FDOWNN*FFLUX* AMUO ) / (PI* ( 1 . -SBARN*RHOS) ) 

( 2 ********************************************************************* 
C * WRITE OUTPUT, ALL MEASURED PARAMETERS WHICH ENTER AS INPUT; AND * 

C *COMPUTED IRRADIANCE AIRR, ABSOLUTE RADIANCE ARAD, AND SURFACE * 

C * REFLECTANCE AT GROUND LEVEL. * 

0 ********************************************************************* 
WRITE ( 6 , 100 6) MTIME, MWAV, MTAU, MTHETO , MTHET, MPHI, MHGHT, AMINT, MINT, 

1 AIRR, ARAD, RHOS 

WRITE ( 5 6 , 1 1 0 7 ) MT IME , MWAV, MTAU, MTHETO , MTHET , MPHI , MHGHT, AMINT, 

1MINT, FINTN, FDOWNN, FTN, SBARN 
GO TO 9987 

0 ********************************************************************* 
C* NOPT IS 2 * 

C* COMPUTE REFLECTANCE MINT (SEE DOCUMENTATION (ABSTRACT) ) * 

0 ********************************************************************* 

C COMPUTE MINT— FINTN+ (FTN*FDOWNN*RHOS ) / ( 1 (SBARN *RHOS) 

C 

56 CONTINUE 
RHOS-MRHO 

MINT— FINTN+ (FTN*FDOWNN*RHOS) / 

1 ( ( 1 . - ( SBARN*RHOS ) ) ) 

Q ******************************************************************** 

C * COMPUTE IRRADIANCE IN WATTS /METER* *2 /UM, AND ABSOLUTE RADIANCE * 

C * IN WATTS/METER* *2/UM/SR. * 

Q ******************************************************************** 

C COMPUTE IRRADIANCE AIRR- (FDOWNN*FFLUX* MUO )/( 1-SBARN . RHOS) 

C COMPUTE ABSOLUTE RADIANCE— (RHOS*FDOWN*FFLUX*AMU0 ) / (PI* ( l-SBAR*RHOS) ) 

C 

AIRR— (FDOWNN * FFLUX* AMUO) / (1 . -SBARN*RHOS) 

ARAD- (RHOS * FDOWNN* FFLUX* AMUO) / (PI* ( 1 . -SBARN*RHOS) ) 

Q ********************************************************************* 

C * WRITE OUTPUT, ALL MEASURED PARAMETERS WHICH ENTER AS INPUT; AND * 

C *COMPUTED IRRADIANCE AIRR, ABSOLUTE RADIANCE ARAD, AND SURFACE * 

C * REFLECTANCE AT GROUND LEVEL. * 

Q ********************************************************************* 

IF ( IOPT .EQ.l .AND. NOPT.EQ.2 ) AMINT- (MINT*FFLUX*AMUO) /PI 
IF ( IOPT . EQ . 2 .AND. NOPT.EQ.2) AMINT-MINT 

IF ( IOPT . EQ . 2 .AND. NOPT.EQ.2) MINT- (MINT*PI) / (FFLUX*AMU0) 

WRITE (6, 1006) MTIME , MWAV, MTAU, MTHETO , MTHET , MPH I , MHGHT, AMINT, MINT, 
1AIRR, ARAD, RHOS 

WRITE ( 5 6 , 1107) MTIME, MWAV, MTAU, MTHETO , MTHET , MPHI , MHGHT , AMINT , 
1MINT, FINTN, FDOWNN, FTN, SBARN 
9987 CONTINUE 
999 CONTINUE 

WRITE (6, 1007) 

STOP 

9 FORMAT (315) 

10 FORMAT (415, 30A1) 

11 FORMAT (/,' DATE ',13 , ' / ' , 13 , ' / ' , 15, ' TIME ZONE',30A1) 
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13 

21 

22 


23 

24 

25 

1000 

1001 

1002 

1003 

1004 

1005 

1006 

1007 

1008 

1107 

1108 
1109 


FORMAT (25X, ' | 27X, ' FLIGHT LEVEL' ,2 8X, ' | ', 13X, ' GROUND LEVEL' ,/, 25XFIF05870 
1,' | ' , 67X, ' | ' , /, 25X, ' |',67X, ' I',/, FIF05880 

1' TIME',2X, 'WAVELENGTH', 2X, ' OPTICAL' ,2X, ' SOLAR', 2X, FIF05890 

1 ' OBSERVATION' ,2X, 'OBSERVATION' ,2X, 'OBSERVATION' , 12X, ' RELATIVE' , 1X,FIF05900 
1' | ' , 2X, 'DOWNWARD' , 5X, 'UPWARD' , /, 18X, ' THICKNESS' , IX, FIF05910 

1' ZENITH' , 3X, ' ZENITH' , 7X, 'AZIMUTH' , FIF05920 

1 6X, 'HEIGHT', 5X, ' RADIANCE', 2X, ' RADIANCE', IX, ' |', 3X, FIF05930 

1 ' IRRADIANCE' , 3X, ' RADIANCE' , 2X, ' REFLECTANCE' , / ,25X, ' | ' , 3X, ' ANGLE' , FIF05940 
14X, 'ANGLE' , 8X, 'ANGLE' ,23X, ' RADIANCE*PI/F0/U0' , 10X, ' W/M**2/UM/SR' FIF05950 

1, /, 6X, 'MICROMETERS' , 8X, ' | ' , 3X, 'DEG' , 6X, 'DEG' , 10X, 'DEG' , 11X, 'KM' , FIF05960 
14X, ' W/M**2/UM/SR' , 10X, ' |',5X,'W/M**2/UM',/,25X,' | ' , 67X, ' | ' , /25X, FIF05970 

1' | ', 67X, ' | ' ) FIF05980 

FORMAT (15, 1X,F6.3,F6.3,4F7.2,F10.5,3F8.4) FIF05990 

FORMAT (/,' JULIAN DATE ' , 13 ,'/',l5,' TIME ZONE', 30A1) FIF06000 

FORMAT ( ' TIME',2X, ' WAVELENGTH', 2X, ' OPTICAL', 2X, ' SOLAR', 2X, FIF06010 

1 ' OBSERVATION' ,2X, 'OBSERVATION' , 2X, ' OBSERVATION' , 12X, ' RELATIVE' , 5X,FIF06020 
l'LSUBO ',2X,'FSUBD ',/,18X, ' THICKNESS', IX, FIF06030 

1' ZENITH' , 3X, ' ZENITH' , 7X, 'AZIMUTH' , FIF06040 

1 6X, 'HEIGHT' , 5X, 'RADIANCE' ,2X, 'RADIANCE' ,2X, FIF06050 

1 ' SUPERPRIME', IX, 'SUPERPRIME', 2X, ' T ' , 2X, ' S ' , / , 2 9X, ' ANGLE' , FIF06060 

14X, 'ANGLE' , 8X, 'ANGLE' , 28X, ' RADIANCE*PI/F0/U0' , /, FIF06070 

16X, 'MICROMETERS' , 13X, 'DEG' , 6X, 'DEG' , 10X, 'DEG' , 11X, 'KM' ,2X, FIF06080 

1 'W/M* *2/UM/SR') FIF06090 

FORMAT (3 (8X,F6.4,F11.4,2X) ) FIFO 6100 

FORMAT (' MEASURED WAVELENGTH (MWAV) ',F6.3,' IS OUT OF RANGE. THE FIFO 6110 
1DATA POINT #' , 14, ' IS NOT PROCESSED. RANGE IS 0.486-2.189') FIF06120 

FORMAT (' MEASURED WAVELENGTH (MWAV) ',F6.3,' IS OUT OF RANGE . THE FIF06130 
1DATA POINT #' , 14, ' IS NOT PROCESSED. RANGE IS 0.486-2.189') FIF06140 

FORMAT ( 6X, 13F6.0) FIF06150 

FORMAT (4X, 1916) FIF06160 

FORMAT (1 IX, F7. 3, 19X, F5 . 2, 22X, F5 . 0, /, 5X, F7 . 4, 11X, F6 . 4) FIF06170 

FORMAT (132A1) FIF06180 

FORMAT (10E12. 4) FIF06190 

FORMAT (3X, 13F6.5) FIF06200 

FORMAT (I5,F9.3,5X,F5.2, IX, ' | ' , 2X, F5 . 2, 3X, F6 . 2 , 7X, F6 . 2 , 7X, FIFO 6210 

1 F7.2,3X,F8.2,4X,F8.4,1X, ' | ' , 3X, F8 . 0, 3X, F7 . 1, 5X, F7 . 3) FIF06220 

FORMAT (25X, ' | ' , 67X, ' | ' , /25X, ' | ' , 67X, ' | ' ) FIF06230 

FORMAT (2F6. 3) FIF06240 

FORMAT (I5,F9.3,5X,F5.2, 6X, F5 . 2, 3X, F6 . 2, 7X,F6.2,7X, FIF06250 

1 F6.2,3X,F8.4,6X,F6.4,4X,F6.4,3X,F6.4,3X,F6.4,3X,F6.4,3X,F6.4) FIFO 6260 
FORMAT ( 45X, ' INTERPOLATED RADIATION PARAMETERS') FIF06270 

FORMAT (43X, ' ' ) FIF06280 

END FIF06290 
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Q it^******************************************************************* 

C * SUBROUTINE READ IN * 

Q ********************************************************************* 

C THIS SUBROUTINE READS THE DATA FROM UNIT NUMBER 5 TO WRITE ON UNIT 
C NUMBER 6. 

C 

SUBROUTINE READ IN 
CHARACTER * 1 LINE (132) 

WRITE (6,200) 

DO 10 1-1,3 
READ (5, 100) LINE 
WRITE (6, 100) LINE 
10 CONTINUE 

READ (5, 101) INUM, (LINE ( I ) , 1-6, 132 ) 

WRITE (6, 101) INUM, (LINE (I) , 1=6, 132) 

READ (5, 100) LINE 
WRITE (6, 100) LINE 
DO 20 1=1, INUM 
READ (5, 100) LINE 
WRITE (6, 100) LINE 
20 CONTINUE 

WRITE (6, 201) 

REWIND 5 

100 FORMAT (132A1) 

101 FORMAT (15, 127A1) 

200 FORMAT (/, IX, 'START OF INPUT DATA AS READ FROM UNIT NUMBER 5',/) 

201 FORMAT (/, IX, ' END OF INPUT DATA SET AS READ FROM UNIT NUMBER 5',/) 
RETURN 

END 

Q **************************★★★*★*★*★*★*★★★*********★*★**★*******★★** 

C * SUBROUTINE FINDW * 

q **★*★**★*★**★*★*★*★*★★**★★★★★*★**★★**★*★★★*★★★★*★*★★★***★******★★★* 
SUBROUTINE FINDW (MWAV, II, INUM, KOPT, IMM,W1,W2) 

REAL* 4 MWAV,WAV(8) / . 486, .587, .639, .663, .837, .845,1.663,2.189/ 

C THIS SUBROUTINE PICKS TWO WAVELENGTHS W1 AND W2 BETWEEN WHICH 
C MEASURED WAVELENGTH MWAV LIES. THE SUBROUTINE PRINTS A ERROR MESSAGE 
C IF THE MWAV DOES NOT FALL IN THE RANGE OF .486 - 2.189 UM.THE 
C CONTROL IS RETURNED TO MAIN PROGRAM WITH VALUE OF LOPT. 

C 

IK = 7 

DO 10 IJ— 1, IK 

IF (MWAV. LT. WAV (1) .OR. MWAV .GT. WAV(IK+1)) GO TO 20 
IF (MWAV.GE.WAV (I J) .AND. MWAV .LE. WAV(IJ+1)) GO TO 30 
10 CONTINUE 
GO TO 20 
30 II-IJ 

IMM = IJ+1 
W1 = WAV ( I J) 

W2 = WAV (IJ+1) 

KOPT-1 

RETURN 

20 WRITE (6, 2) MWAV, INUM 
KOPT-0 

2 FORMAT (' MEASURED WAVELENGTH (MWAV) ',F6.3,' IS OUT OF RANGE. THE 
1DATA POINT #' , 14, ' IS NOT PROCESSED. ACTUAL RANGE .486-2.189 UM' ) 
RETURN 
END 
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Q A******************************************************************** 

C * SUBROUTINE INTSFX * 

Q ********************************************************************* 

SUBROUTINE INTSFX (RR) 

THIS SUBROUTINE INTERPOLATES VALUES OF RR FOR 365 DAYS OF YEAR. 

RR IS SUN EARTH DISTANCE. 

REAL * 4 RR (365) , DAY (365) 

REAL * 4 AMON(25) /l. , 15. , 32 . , 46. , 60 . , 74 ., 91 . , 106 . , 

1121.. 135.. 152.. 166.. 182.. 196.. 213.. 227., 

1 224., 258., 274., 288., 305., 319., 325., 349., 365./ 

REAL * 4 R(25) /. 9832, .9836, .9853, .9878, 

1.9909. . 9945, .9993,1.0033,1.0076, 

1 1.0109,1.0140,1.0158,1.0167,1.0165,1.0149,1.0128,1.0092,1.0057, 

1 1.0011, .9972, .9925, .9892, .9860, .9843, .9843/ 

DO 100 1=1,365 
DAY (I) =1 
100 CONTINUE 
IJK-0 
NUM=0 

CALL INTERP (365,25, DAY, AMON , R, RR, I JK, NUM, LOPT) 

RETURN 
END 

********************************************************************* 

* SUBROUTINE INTHGH * 

********************************************************************* 

SUBROUTINE INTHGH (LTAU, LTHE, LPHI, LTHET0 , INTW, MHGHT, PITW, INTWH, 

1 PITWH) 

THIS SUBROUTINE INTERPOLATES THE VALUES OF INTW AND PITW FOR THE 
MEASURED HEIGHT MHGHT 

REAL * 4 INTW(2, 9, 4, 3, 13, 19) , INTWH(2, 9, 4, 13, 19) ,PITW(2, 13, 4, 3) , 

1 PITWH (2 , 13, 4 ) , HGHT ( 3 ) /80., 4.5, .45/, MHGHT 

REAL * 4 INTWH1 (2, 9,4, 13, 19) ,PITWH1 (2, 13, 4) 

REAL * 4 INTWH2 (2, 9, 4, 13, 19) ,PITWH2 (2, 13, 4) 

C 

C INTERPOLATES IF MEASURED HEIGHT IS LESS THAN OR EQUAL TO .45KM 
C (SEE DOCUMENTATION EQU. 21A, 2 IB) 

C 

IF ( MHGHT .LE. HGHT (3)) THEN 
DO 100 I WAV =1,2 
DO 101 ITAU = 1 , LTAU 
DO 102 ITHE = 1, LTHE 

PITWH (IWAV, ITHE, ITAU) =( (PITW ( IWAV, ITHE, ITAU, 3) * MHGHT) + 

1 (HGHT (3) -MHGHT) ) /HGHT (3) 

DO 103 IPHI = 1 , LPHI 

DO 104 ITHET0 = 1, LTHET0 

INTWH ( IWAV, ITHET0 , ITAU, ITHE, IPHI ) = 

1 (INTW (IWAV, ITHET0, ITAU, 3, ITHE, IPHI) * MHGHT) / HGHT (3) 

IF (INTWH (IWAV, ITHET0, ITAU, ITHE, IPHI) . LE . 0.0) 

1 INTWH (IWAV, ITHET0, ITAU, ITHE, IPHI) =1 .E-6 
104 CONTINUE 
103 CONTINUE 
102 CONTINUE 
101 CONTINUE 
100 CONTINUE 
C 
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C INTERPOLATES IF MEASURED HEIGHT LIES BETWEEN .45 AND 4 . 5 KM 
C (SEE DOCUMENTATION 3.6) 

C 

ELSE 

IF ( MHGHT .GT. HGHT ( 3 ) .AND. MHGHT .LT. HGHT (2) ) THEN 


EXTRAPOLATES FOR THE HEIGHT LE . TO .45 KM. 


DO 200 I WAV - 1,2 
DO 201 ITAU - 1 , LTAU 
DO 202 ITHE - 1, LTHE 

PITWH1 (IWAV, ITHE, ITAU) - ( (PITW (IWAV, ITHE, ITAU, 3) * MHGHT) + 

1 (HGHT (3) -MHGHT) ) /HGHT (3) 

DO 203 IPHI - 1, LPHI 
DO 204 ITHET0 - 1, LTHET0 
INTWH1 (IWAV, ITHET0, ITAU, ITHE, IPHI) - 
1 (INTW( IWAV, ITHET0, ITAU, 3, ITHE, IPHI) * MHGHT) / HGHT (3) 

IF (INTWH1 (IWAV, ITHET0, ITAU, ITHE, IPHI) .LE. 0.0) 

1 INTWH1 (IWAV, ITHET0, ITAU, ITHE, IPHI) =1 .E-6 
CONTINUE 
203 CONTINUE 
202 CONTINUE 
201 CONTINUE 
CONTINUE 


EXTRAPOLATES BETWEEN THE HEIGHTS 4 . 5 KM AND 80.0 KM. 

2 - (1. - (EXP (-MHGHT/ 9.) ) ) 

21- (1. - (EXP (-HGHT (2) /9 . ) ) ) 

22- (1. - (EXP (-HGHT (l)/9.))) 

DO 300 IWAV - 1,2 

DO 301 ITAU - 1 , LTAU 
DO 302 ITHE - 1 , LTHE 
PITWH2 (IWAV, ITHE, ITAU) - 
1PITW( IWAV, ITHE, ITAU, 2) + 

1 (( (PITW (IWAV, ITHE, ITAU, 1) - 

1 PITW (IWAV, ITHE, ITAU, 2) ) / (22 -21)) * (2 -21)) 

DO 303 IPHI - 1, LPHI 
DO 304 ITHET0 - 1, LTHET0 
INTWH2 (IWAV, ITHET0, ITAU, ITHE, IPHI) - 
1 INTW ( IWAV, ITHET0, ITAU, 2, ITHE, IPHI) + 

1 (( (INTW (IWAV, ITHET0, ITAU, 1, ITHE, IPHI) - 

1 INTW ( IWAV, ITHET0, ITAU, 2, ITHE, IPHI) ) / (22 -21)) * (2 -21)) 

04 CONTINUE 
03 CONTINUE 
302 CONTINUE 
301 CONTINUE 
CONTINUE 


CHOSSES THE VALUES THAT SHOW NOMINAL ATMOSPHERIC EFFECT. LOWER VALUES 
OF INTENSITY AND HIGHER VALUES OF PIT ARE PICKED. 
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DO 400 IWAV - 1,2 
DO 401 ITAU - 1 , LTAU 
DO 402 ITHE - 1 , LTHE 

IF (PITWHl (IWAV, ITHE, ITAU) . LE . PITWH2 ( IWAV, ITHE, ITAU) ) GO TO 405 
PITWH (IWAV, ITHE, ITAU) - PITWHl ( IWAV, ITHE, ITAU) 

DO 403 IPHI - 1 , LPHI 

DO 404 ITHET0 =* 1, LTHET0 

INTWH (IWAV, ITHET0, ITAU, ITHE, IPHI) = 

1INTWH1 (IWAV, ITHET0, ITAU, ITHE, IPHI) 

404 CONTINUE 
403 CONTINUE 

GO TO 402 

405 CONTINUE 

PITWH (IWAV, ITHE, ITAU) - PITWH2 (IWAV, ITHE, ITAU) 

DO 407 IPHI - 1 , LPHI 

DO 406 ITHET0 - 1, LTHET0 

INTWH (IWAV, ITHET0, ITAU, ITHE, IPHI) = 

1INTWH2 (IWAV, ITHET0, ITAU, ITHE, IPHI) 

406 CONTINUE 

407 CONTINUE 
402 CONTINUE 
401 CONTINUE 
400 CONTINUE 

INTERPOLATE FOR THE HEIGHT BETWEEN 4.5 KM. AND 80.0 KM. 

ELSE 

IF ( MHGHT .GE. HGHT (2) ) THEN 
Z - (1. - (EXP < -MHGHT/ 9.) ) ) 

Zl- (1. - (EXP (-HGHT (2) /9 . ) ) ) 

Z2- (1. - (EXP (-HGHT (l)/9.))) 

C WRITE(6,10001)Z,Z1,Z2 

DO 600 IWAV - 1,2 
DO 601 ITAU - 1 , LTAU 
DO 602 ITHE = 1 , LTHE 
PITWH (IWAV, ITHE, ITAU) « 

IP ITW( IWAV, ITHE, ITAU, 2) + 

1 (( (PITW(IWAV, ITHE, ITAU, 1) - 

1 PITW(IWAV, ITHE, ITAU, 2) ) / (Z2 - Zl) ) * (Z -Zl) ) 

DO 603 IPHI = 1, LPHI 
DO 604 ITHET0 “ 1, LTHET0 
INTWH (IWAV, ITHETO, ITAU, ITHE, IPHI) - 
1 INTW(IWAV, ITHETO, ITAU, 2, ITHE, IPHI) + 

1 (( <INTW(IWAV, ITHETO, ITAU, 1, ITHE, IPHI) - 

1 INTW (IWAV, ITHETO, ITAU, 2, ITHE, IPHI) ) / (Z2 - Zl) ) * (Z -Zl) ) 

604 CONTINUE 
603 CONTINUE 
602 CONTINUE 
601 CONTINUE 
600 CONTINUE 
ELSE 

END IF 
END IF 
END IF 
RETURN 
END 
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c 

c 

c 

c 

c 

c 

c 


********************************************************************* 

* SUBROUTINE INTERP * 

********************************************************************* 

SUBROUTINE INTERP (I,M, XI, X, Y, Yl r IJK, INUM, LOPT) 

THIS SUBROUTINE IS A GENERAL PURPOSE ROUTINE AND INTERPOLATES 
LINEARLY. VALUE OF Y1 IS INTERPOLATED FOR XI. NO EXTRAPOLATION IS 
ALLOWED . 

REAL *4 XI (I) , X (M) , Y (M) , Y1 (I) 

DO 290 LK = 1,1 
XBAR=X1 (LK) 

LL-M-1 

DO 230 IL«1,LL 

IF (XBAR . LT . X ( 1 ) ) GO TO 300 

IF (XBAR . GT .X (LL+1) ) GO TO 300 

IF (XBAR .GE.X(IL) . AND . XBAR . LE . X(IL+1)) GO TO 250 
GO TO 230 
250 PPHI-X(IL) 

SPHI«X(IL+1) 

PINTEN=Y(IL) 

SINTEN=Y ( IL+1 ) 

30 Y1 (LK) “PINTEN+ ( (SINTEN-PINTEN) * ( (XBAR-PPHI) / (SPHI-PPHI) ) ) 

GO TO 290 
230 CONTINUE 
290 CONTINUE 
LOPT=l 
RETURN 

THE PROGRAM DOES NOT ALLOW ANY EXTRAPOLATION .IT WILL PRINT AN ERROR 
MESSAGE .IJK IDENTIFIES THE ARRAY BEING SENT FOR INTERPOLATION. 


C 
C 
300 


LOPT-O 
IF (IJK 
IF (IJK 
IF (IJK 
IF (IJK 
IF (IJK 
IF (IJK 
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. EQ . 1 ) WRITE (6, 1) XBAR, INUM 
. EQ . 2 ) WRITE (6, 1) XBAR, INUM 
. EQ . 3 ) WRITE (6, 2) XBAR, INUM 
. EQ . 4 ) WRITE (6, 3) XBAR, INUM 
. EQ . 5 ) WRITE (6, 3) XBAR, INUM 
. EQ . 6 ) WRITE (6, 4) XBAR, INUM 

1 FORMAT (IX, ' MEASURED SOLAR ZENITH ANGLE (MTHET0) ', F5 . 1, ' IS OUT OFFIF08900 

1 RANGE. ACTUAL RANGE 10 -78 DEGREES. THE DATA POINT #',I4,' IS NOT FIF08910 
1 PROCESSED.') FIF08920 

2 FORMAT ( IX, ' MEASURED OBSERVATION AZIMUTH ANGLE (MPHI) ', F6 . 1, ' IS OUFIF08930 

IT OF RANGE. ACTUAL RANGE 0 -180 DEGREES. THE DATA POINT #' , 14, ' IS NFIF08940 
10T PROCESSED.') FIF08950 

3 FORMAT ( IX, ' MEASURED OBSERVATION ZENITH ANGLE (MTHETA) ', F5 . 1 , ' IS OFIF08960 

1UT OF RANGE. ACTUAL RANGE 0 -72 DEGREES. THE DATA POINT #',I4,' IS NFIF08970 
10T PROCESSED.') FIF08980 

4 FORMAT (IX,' MEASURED OPTICAL THICKNESS (MTAU) ',F5.3,' IS OUT OF RFIF08990 

1ANGE . ACTUAL RANGE 0. 0-1.0 . THE DATA POINT #' , 14, ' IS NOT PROCESFIF0 9000 
1SED . ' ) FIF09010 

RETURN FIF09020 

END FIF09030 

SUBROUTINE INTEXP (X, Y, MTAU, YY, LTAU, IGNORE) FIF09040 

C FIF09050 

C THIS SUBROUTINE IS A GENERAL PURPOSE ROUTINE AND INTERPOLATES FIF09060 

C EXPONANTIALLY. IGNORE SENDS VALUE OF 0 IF SUBROUTINE INTEXP IS FIF09070 

C USED AND 1 IF IT FINDS THE FUNCTION TO BE LINEAR. FIF09080 

C FIF09090 

REAL * 4 X (LTAU) ,Y (LTAU) ,MTAU,YY(1) , TOL, C FIF09100 

COMMON/FF/XS (3) , YS (3) FIF09110 

EXTERNAL FMIN FIFO 9120 

C FIF09130 
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IGNORE=0 
IF (MTAU.LE. 
MM=1 
MN=3 


0.25) THEN 


10 


IF MEASURED OPTICAL THICKNESS IS GREATER THAN 0.25 THEN FUNCTIONS 
OF OPTICAL THICKNESSES OF 0.25, 0.50 AND 1.00 ARE USED. 

ELSE 
MM=2 
MN-4 
END IF 
IS=0 

DO 10 IJ - MM, MN 

IS-IS+1 

XS (IS)-X(IJ) 

YS(IS)-Y(IJ) 

CONTINUE 


SYSTEMS SUBROUTINE 
FUNCTION FMIN . 


IMSL (ZXGSN) IS CALLED TO GET THE MINIMAL OF 


C 

C 

C 

C 

C 


VALUE OF IGNORE IS SET TO 0 AND IF MEASURED OPTICAL THICKNESS IS LESS FIF09140 
THEN OR EQUAL TO 0.25 THEN THE FUNCTIONS OF OPTICAL THICKNESSES OF FIF09150 
0.0 ,0.25, 0.50 ARE USED. FIF09160 
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POINT OF THE INTERVAL IN WHICH MINIMUM OF FMIN IS TO BEFIF09400 

FIF09410 
FIF09420 
FIF09430 
FIF09440 
FIF09450 
FIF09460 
FIF09470 
FIF09480 
FIF09490 
FIF09500 
FIF09510 
FIF09520 
FIF09530 
FIF09540 
FIF09550 
FIF09560 
FIF09570 
FIF09580 
FIF09590 
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IF FUNCTION IS NOT LINEAR , EXPONANTIAL INTERPOLATION IS USED, AND INTEFIF09630 
INTERPOLATED VALUE FOR MTAU IS YY(1) IS COMPUTED FOR THE FUNCTION. FIFO 9640 

FIFO 9650 

BB= (YS (1) -YS (2) ) / ( (EXP (~C*XS (1) ) ) - (EXP (-C*XS (2) ) ) ) FIF09660 


A IS LOWER END- 
LOCATED . 

B IS THE UPPER END-POINT OF THE INTERVAL. 

TOL IS LENGTH OF THE FINAL SUBINTERVAL CONTAINING THE MINIMUM. 

C IS OUTPUT , THE APPROXIMATE MINIMUM OF THE FUNCTION FMIN ON THE 
ORGINAL INTERVAL (A, B) 

IER RETURNS THE VALUE OF 0 IF NO ERROR IS FOUND . 


A - 
AA ■ 
B - 
TOL ! 


0.01 
■ 0.01 
10 

■1 .E-4 


CALL ZXGSN (FMIN, A, B, TOL, C, IER) 

THE IF STATEMENT BELOW DECIDES IF THE FUNCTION IS LINEAR OR NOT 
IF THE FUNCTION IS LINEAR THE VALUE OF IGNORE IS SET TO 1 AND 
CONTROL IS SEND BACK TO MAIN PROGRAM. 


IF ( (C-AA) 
IGNORE=l 
RETURN 
ELSE 


.LT. .001)THEN 


AA= Y S ( 1 ) - ( BB * (EXP (-C*XS (1) ) ) ) 
YY ( 1 ) =AA+ (BB * (EXP (-C*MTAU) ) ) 
END IF 
RETURN 
END 
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FUNCTION FMIN(C) 


FIF09720 

FIF09730 

FUNCTION FMIN IS FUNCTION TO BE MINIMIZED. 

FIF09740 

FIF09750 

REAL C 


FIF09760 

COMMON/FF/XS (3) , YS (3) 


FIF09770 

P = (YS (1) - YS (2) ) / 

(YS (3) - YS (2) ) 

FIF09780 

PF - EXP ( -C * XS(1) ) - 

EXP (-C*XS (2) ) 

FIF09790 

PG - EXP ( -C * XS (3) ) - 

EXP (-C*XS (2) ) 

FIF09800 

PP-PF/PG 


FIF09810 

FMIN-(P-PP) **2 


FIF09820 

RETURN 


FIF09830 

END 


FIF09840 
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